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1.0  ABSTRACT 


The  method  described  here  is  a  higher-order  three-dimensional  surface 
panel  method.  It  utilizes  four-sided  panels  on  the  actual  surface  of  the 
body  about  which  flow  is  to  be  computed  and  is  thus  applicable  to  arbitrary 
configurations.  In  contrast  to  a  first-order  panel  method  that  uses  flat 
panels,  usually  with  constant  source  and/or  vorticity  density,  the  present 
method  uses  curved  panels  of  second  degree  with  linearly  varying  source  and 
vorticity  density.  The  mathematical  problem  solved  is  that  of  incompressible 
potential  flow,  but  the  calculation  can  be  extended  to  subsonic  flow  by  means 
of  well-known  compressibility  corrections. 

This  report  presents  a  discussion  of  the  techniques  adopted  and  some 
alternatives.  A  rigorous  derivation  of  the  higher-order  formulation  is 
carried  out.  All  formulas  and  logic  required  by  the  method  are  included, 
except  that  the  report  on  the  corresponding  first-order  method  is  relied  on 
for  formulas  and  procedures  common  to  the  two  techniques.  First-order  and 
higher-order  calculated  results  are  compared  with  each  other  and  with  exact 
solutions. 
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3.0  INTRODUCTION 


Recent  years  have  witnessed  the  acceptance  of  so-called  "panel  methods" 
for  calculating  three-dimensional  potential  flows  to  such  an  extent  that  they 
are  now  considered  to  be  routine  tools  of  aerodynamic  and  hydrodynamic  design. 
While  some  supersonic  work  has  been  done,  by  far  the  chief  emphasis  in  panel 
method  development  has  been  directed  to  the  solution  of  Laplace's  equation 
for  incompressible  flow.  The  results  obtained  are  generalized  to  subsonic 
flows  by  the  application  of  certain  compressibility  corrections.  Panel  methods 
may  be  grouped  into  two  broad  categories.  Surface  panel  methods  place  panels 
and  satisfy  normal  velocity  boundary  conditions  on  the  surface  of  the  body 
about  which  flow  is  to  be  computed.  (Actually  the  numerical  surface  is  a 
close  approximation  to  the  body  surface,  which  it  approaches  in  the  limit  of 
infinite  panel  number.)  This  type  of  panel  method  is  unrestricted,  and  exact 
in  the  sense  that  it  is  applicable  to  completely  arbitrary  bodies  and  that 
incompressible  flow  may  be  calculated  to  any  degree  of  precision  by  increasing 
the  panel  number.  On  the  contrary,  mean-surface  panel  methods  place  panels 
interior  to  the  body  in  question,  e.g.,  on  the  axis  of  a  body  of  revolution 
or  on  the  mean  surface  of  a  wing.  This  type  of  panel  method  is  approximate 
and  is  applicable  to  slender  bodies,  and  to  cases  where  the  flow  perturbation 
velocities  are  small.  Both  types  of  method  have  their  uses.  Mean-surface 
panel  methods  are  used  for  preliminary  design,  while  surface  panel  methods 
are  used  for  final  detailed  calculations.  It  is  with  incompressible  surface 
panel  methods  that  the  present  report  is  concerned. 

Panel  methods  arise  as  numerical  discretizations  of  the  theoretical 
technique  that  may  be  denoted  the  surface  singularity  approach.  In  this 
procedure  the  surface  of  the  body  about  which  flow  is  to  be  computed  is 
imagined  to  be  covered  with  distributions  of  flow  singularity:  source,  dipole, 
vorticity,  and  the  variable  strengths  of  these  distributions  are  adjusted  to 
satisfy  the  normal  velocity  boundary  condition  on  the  body  surface  and  any 
auxiliary  conditions,  such  as  Kutta  conditions  along  a  wing  trailing  edge. 

In  particular  the  normal  velocity  boundary  condition  leads  to  an  integral 
equation  in  the  singularity  strengths.  In  the  numerical  implementation  of 
this  technique  the  first  step  must  be  to  discretize  the  body  surface.  In 
three-dimensional  problems  the  body  surface  is  represented  by  a  large  number 


of  small  four-sided  surface  elements  or  "panels"  which  is  the  origin  of  the 
name  "panel  methods." 

There  are  various  possibilities  for  distributing  types  of  singularities 
on  the  surface  panels  and  different  investigators  have  made  different  choices 
for  their  methods.  However,  it  is  conmon  to  all  methods  that  one  particular 
type  of  singularity  is  assigned  primary  responsibility  for  satisfying  the 
normal  velocity  boundary  condition  and  thus  may  be  called  the  basic  singular¬ 
ity.  Other  types  of  singularity  have  primary  responsibility  for  auxiliary 
conditions  and  may  be  called  auxiliary  singularities.  (Clearly  all  conditions 
must  ultimately  be  satisfied  simultaneously  but  the  above  is  a  valid  distinc¬ 
tion.)  It  has  become  customary  to  categorize  a  panel  method  according  to  its 
basic  singularity.  Thus,  reference  is  made  to  source  methods,  dipole  methods, 
and  vortex  methods.  This  report  is  concerned  with  generalizations  of  the 
Hess  panel  method  (refs.  1  and  2),  which  is  in  such  common  use  throughout  the 
country  that  it  may  fairly  be  called  the  standard  existing  procedure.  The 
Hess  technique  is  a  source  method,  and  the  descriptions  below  apply  to  that 
type  of  method,  although  parallel  descriptions  for  the  other  methods  would 
differ  only  in  detail. 

The  panel-method  representation  of  a  three-dimensional  body  is  shown  in 
figure  1.  On  each  panel  a  control  point  is  selected  where  the  boundary  condi¬ 
tion  of  zero  normal  velocity  is  to  be  applied  and  where  surface  velocities  are 
ultimately  calculated.  A  "matrix  of  influence  coefficients"  is  then  calculated. 
This  consists  of  the  complete  set  of  velocities  induced  by  the  panels  at  each 
others'  control  points  for  certain  nominal  variations  of  the  singularity 
strengths,  e.g.,  constant  at  unit  value,  linearly  varying  with  unit  deriva¬ 
tive,  etc.  On  nonlifting  portions  of  the  configuration  (figure  1)  only  source 
singularity  is  employed.  On  lifting  portions  both  source  and  either  dipole 
or  vorticity  singularities  are  used.  It  can  be  shown  (ref.  1)  that  surface 
distributions  of  dipole  and  of  vorticity  are  equivalent.  The  dipole  has 
been  used  most  commonly  in  three  dimensions  because  it  is  a  scalar  rather  than 
a  vector  and  because  it  automatically  satisfies  the  appropriate  conservation 
theorems  that  must  be  applied  to  vorticity  distributions  (essentially,  zero 
divergence  of  vorticity).  From  the  basic  standpoint  a  lifting  portion  of  the 
configuration  is  defined  as  one  having  a  sharp  trailing  edge  where  a  Kutta 
condition  is  to  be  applied  and  from  which  issues  a  wake  of  trailing  vorticity 
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(figure  1).  Nonlifting  portions  are  characterized  by  the  absence  of  a  trail¬ 
ing  edge  and  its  associated  vortex  wak?-.  The  prototypes  of  lifting  and  non¬ 
lifting  portions  are  a  wing  and  a  fuselage,  respectively.  On  the  wake  only 
dipole  (vorticity)  singularity  is  needed.  However,  from  the  panel-method 
standpoint  there  is  a  change  of  emphasis.  A  lifting  portion  is  characterized 
by  one  having  bound  vorticity  whose  strength  is  adjusted  to  satisfy  certain 
conditions,  while  a  nonlifting  portion  has  no  vorticity  and  no  conditions. 

A  panel  method  approximates  the  integral  equation  that  expresses  the  zero 
normal -velocity  boundary  condition  by  a  set  of  linear  algebraic  equations  for 
the  values  of  singularity  strengths  on  the  panels.  In  lifting  cases  the  Kutta 
condition  along  the  trailing  edge  (figure  1)  yields  a  small  number  of  addi¬ 
tional  equations.  On  each  panel  one  point,  which  is  designated  the  control 
point,  is  selected,  and  the  normal  velocity  is  required  to  be  zero  there. 

Each  panel  has  an  independent  value  of  source  strength.  This  provides  a  num¬ 
ber  of  unknown  parameters  equal  to  the  number  of  panels  and,  more  particularly 
equal  to  the  number  of  linear  equations  expressing  the  condition  of  zero  normal 
velocity.  The  Kutta  condition  is  discretized  by  applying  it  at  a  finite  number 
of  locations  along  the  trailing  edge.  The  dipole  (vorticity)  distribution 
over  the  surface  is  then  specialized  to  be  expressible  in  terms  of  a  number 
of  adjustable  parameters  equal  to  the  number  of  locations  at  which  the  Kutta 
condition  is  applied.  Once  this  has  been  accomplished,  there  are  the  same 
number  of  unknown  values  of  singularity  as  there  are  conditions  to  be  satisfied. 
In  nonlifting  flows,  no  Kutta  condition  is  applied,  and  the  source  strengths 
are  calculated  to  yield  zero  normal  velocity  at  the  control  points.  This  sum¬ 
mary  is  necessarily  brief.  Details  are  contained  in  references  1  and  2. 

In  the  standard  form  of  the  Hess  panel  method  in  both  its  nonlifting 
(ref.  2)  and  lifting  (ref.  1)  versions,'  the  surface  panels  of  figure  1 
are  plane  quadrilaterals  and  the  source  strength  on  each  panel  is  assumed 
constant  —  a  step  function  from  element  to  element.  (The  dipole  variation, 
however,  is  more  elaborate,  (ref.  1).)  This  standard,  which  is  referred  to 
as  the  base  method,  has  proved  satisfactory  in  numerous  applications  by 
investigators  in  industry,  government,  and  academia.  However,  experience  has 
indicated  that  certain  improvements  would  be  desirable.  Probably  the  princi¬ 
pal  disadvantage  of  the  Hess  method  or  any  three-dimensional  panel  method  is 
that  each  case  is  relatively  expensive.  This  expense  consists  of  two  parts: 


the  human  labor  necessary  to  prepare  the  rather  large  amount  of  input  data 
and  the  computer  cost  required  to  carry  out  the  calculations.  Ironically, 
the  problem  of  computing  time  has  been  intensified  by  the  very  success  of  the 
method  in  application.  Namely,  its  success  has  led  to  its  use  in  more  and 
more  challenging  problems.  The  result  is  that  some  of  the  most  frequent  users 
of  the  method  have  found  themselves  employing  ever  larger  numbers  of  surface 
panels  to  represent  complicated  configurations  accurately.  This  trend  has  been 
in  evidence  for  some  time.  For  a  long  while,  the  largest  panel  number  used 
was  1000.  However,  recent  experience  has  greatly  increased  this  number.  The 
largest  case  run  to  date  employed  3500  surface  panels  on  each  half  of  the 
configuration  (symmetry  was  employed  so  that  the  order  of  the  influence  matrix 
was  3500).  Computing  time  for  a  panel  method  increases  as  the  square  of  the 
panel  number  or  faster.  (If  a  direct  solution  of  the  linear  equations  is  used, 
the  time  can  vary  as  the  cube  of  the  panel  number.)  Thus  large  panel  numbers 
require  very  long  computation  times  and  thus  very  large  costs.  Moreover,  the 
accuracies  obtained  are  not  always  what  is  desired.  Accordingly,  there  is  a 
need  for  an  increase  in  speed  of  computation  with  no  decrease  in  accuracy,  or 
better  still,  with  an  increase  in  accuracy.  This  has  been  the  aim  of  the  work 
described  in  this  report. 

To  reduce  the  other  major  expense,  i.e.,  the  man-hours  required  for 
inputting  a  body  to  a  panel  method,  a  so-called  "geometry  package"  has  been 
constructed  (refs.  3  and  4).  This  routine,  which  is  a  major  program  in  its 
own  right,  enables  a  user  of  the  method  to  input  a  relatively  small  number 
of  points  defining  his  particular  body.  Moreover,  relatively  little  care 
will  be  required  in  selecting  the  particular  points  to  input.  The  "geometry 
package"  will  then  enrich  the  point  number,  redistribute  the  points  as  needed 
according  to  one  of  a  number  of  algorithms,  compute  intersections  of  various 
portions  of  the  configuration,  etc.  For  certain  frequently-occurring  cases 
the  procedure  will  be  completely  automatic.  For  others,  there  will  still  be 
a  "man  in  the  loop,"  but  his  labor  will  be  greatly  reduced.  Figure  2a  shows 
a  very  sparse  representation  of  a  wing- fusel  age  that  was  input  to  the  geometry 
package  which  then  produced  the  enriched  panelling  of  figure  2b.  Included  in 
the  process  is  the  calculation  of  the  wing- fusel  age  intersection  curve.  This 
geometry  package  has  been  combined  with  the  higher-order  program  of  the  present 
report. 
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Both  logically  and  mathematically,  the  present  method  is  an  extension  of 
the  base  method  of  reference  1.  Indeed  many  formulas  and  procedures  are  common 
to  the  two  programs.  If  the  present  report  were  to  be  made  self-contained, 
it  would  have  to  contain  a  very  large  percentage  of  the  material  presented  in 
reference  1.  This  would  not  only  be  repetitious  but  would  make  the  present 
report  extremely  lengthy.  Accordingly,  it  is  assumed  that  a  reader  inter¬ 
ested  in  details  has  reference  1  available,  and  reference  is  freely  made  to 
equations  and  sections  of  that  report.  However,  it  is  possible  to  obtain  a 
general  familiarity  with  the  present  method  without  using  reference  1. 


4.0  GENERAL  DESCRIPTION  OF  THE  HIGHER-ORDER  PANEL  METHOD 


The  purpose  of  the  work  described  in  this  report  is  to  increase  the  speed 
and  accuracy  of  the  Hess  panel  method  by  use  of  a  higher-order  formulation. 

The  higher-order  formulation  replaces  the  standard  base  method's  assumption 
of  plane  surface  panels  having  constant  values  of  source  strength  (and  a  quadratic 
dipole  strength)  by  an  assumption  of  curved  panels  with  variable  source  and 
dipole  strength.  A  direct  analogy  of  this  approach  has  already  been  applied 
with  great  success  to  the  two-dimensional  case  (ref.  5)  and  the  axisymmetric 
case  (ref.  6).  The  concept  of  the  higher-order  approach  was  first  put  for¬ 
ward  in  reference  5,  where  it  was  shown  that  the  effect  of  surface  curvature 
and  the  effect  of  the  variation  of  surface  singularity  strength  are  of  the 
same  order  of  magnitude.  Thus,  inclusion  of  one  without  the  other,  as  is 
done  for  instance  by  several  investigators  in  two  dimensions,  is  inconsistent 
and,  in  general,  does  not  lead  to  an  improvement  in  accuracy  over  the  base 
method.  Use  of  the  consistent  higher-order  formulation  has  yielded  large 
increases  in  both  computational  speed  and  accuracy  for  two-dimensional  and 
axisymmetric  flows  (refs.  5,  6  and  7).  The  present  three-dimensional  higher- 
order  formulation  yields  comparable  gains  in  speed  and  accuracy. 

Conceptually,  the  essence  of  constructing  a  higher-order  panel  method 
consists  of  generating  an  altered  "matrix  of  influence  coefficients"  (previous 
section).  All  other  aspects  of  the  program,  e.g.,  input  of  the  body  geometry, 
solution  of  the  linear  equations,  calculation  of  surface  velocities  and  pres¬ 
sures,  remain  the  same.  Specifically,  the  formulas  that  in  the  base  method 
give  the  velocity  induced  at  an  arbitrary  control  point  by  a  flat  panel  with 
a  constant  source  strength  must  be  replaced  in  a  higher-order  method  by 
formulas  that  give  the  velocity  induced  at  an  arbitrary  control  point  by  a 
curved  panel  with  a  variable  source  strength.  Similarly  the  formulas  that  in 
the  base  method  give  the  velocity  due  to  a  quadratic  dipole  strength  over  a 
flat  panel  are  replaced  by  those  appropriate  to  a  curved  panel.  Determination 
of  the  form  of  the  panel  shape  and  the  variation  of  the  singularity  strength 
has  been  part  of  the  problem.  It  turned  out  that  derivation  of  new  velocity 
formulas  was  only  half  the  required  task.  A  considerable  amount  of  new  logic 
was  also  required. 
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The  starting  point  of  the  present  work  is  the  consistency  analysis  of 
reference  8,  which  considers  a  surface-source  distribution.  It  is  shown, 
therein,  that  consistent  approximations  in  ascending  order  of  accuracy  are: 

(1)  flat-panel  constant- source  (base  method),  (2)  paraboloidal-panel  linear- 
source,  and  (3)  cubic-panel  quadratic-source.  In  other  words,  a  consistent 
approach  always  uses  a  source  polynomial  one  degree  less  than  the  panel  poly¬ 
nomial.  Reference  8  also  gives  the  form  of  the  expansion  and  the  number  of 
basic  integrals  over  the  panel  that  are  required  by  each  order  of  approxima¬ 
tion.  Thus,  the  flat-panel  constant-source  approach  requires  a  single  integral. 
The  paraboloidal-panel  linear-source  approach  requires  five  additional  inte¬ 
grals  -  a  total  of  six.  Finally,  the  cubic-panel  quadratic-source  approach 
requires  a  total  of  23  integrals.  This  last  is  felt  to  be  too  complicated  to 
be  feasible  and  thus  the  "higher-order"  procedure  adopted  is  the  paraboloidal- 
panel  linear-source.  This  is  the  direct  analogy  of  the  formulation  that  proved 
so  effective  in  the  two-dimensional  and  axi symmetric  cases  (refs.  5  and  6).  For 
completeness  the  expansions  derived  in  reference  8  are  presented  in  section  5.1, 
and  thus  the  present  report  supercedes  reference  8. 

A  second  consistency  analysis  is  necessary  for  the  dipole  or  vorticity 
distribution  required  to  produce  lift.  From  specific  two-dimensional  formulas, 
as  well  as  general  three-dimensional  considerations,  it  is  expected  that  a 
consistent  vorticity  distribution  has  the  same  degree  as  the  source  distribu¬ 
tion.  Thus,  according  to  the  relation  between  dipole  and  vorticity  distribu¬ 
tions  (ref.  1),  the  consistent  dipole  distribution  has  a  degree  one  higher 
than  that  of  the  source  distribution.  However,  when  the  analysis  is  performed 
for  the  integrated  effect  of  a  single  panel,  it  turns  out  that  the  quadratic 
dipole  terms  are  of  the  same  order  as  cubic  geometry  terms  (as  is  the  case  for 
the  source  distribution).  This  unfortunate  behavior  is  due  to  the  so-called 
"edge  vortex."  As  derived  in  reference  1,  a  variable  dipole  distribution  over 
a  finite  area  is  equivalent  to  the  sum  of  a  vorticity  distribution  over  the 
area  with  strength  equal  to  the  gradient  of  the  dipole  strength  plus  a  concen¬ 
trated  vortex  filament  around  the  edge  of  the  area  with  variable  strength 
equal  to  the  local  dipole  strength.  It  is  due  to  the  edge  vortex  that  the 
cubic  geometry  terms  enter  to  a  lower  order  than  expected.  Clearly  the  actual 
wing  of  figure  1  does  not  have  concentrated  vortices  along  panel  edges.  Their 
presence  is  due  to  the  fact  that  formulas  are  derived  for  one  isolated  panel. 

The  effects  of  the  edge  vortex  on  a  panel  are  presumably  cancelled  by  equal 
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and  opposite  edge-vortices  on  adjacent  panels.  If  the  geometry  and  the  dipole 
strength  are  continuous  across  panel  edges,  this  cancellation  is  exact. 

Otherwise  it  is  approximate.  However,  the  edge-vortex  effects  enter  the  basic 
expansion  formulas  in  a  way  that  makes  it  difficult  to  cancel  them  analytically. 

The  procedure  adopted  to  avoid  the  above  difficulty  is  the  use  of  a 
vorticity  distribution,  rather  than  a  dipole  distribution,  on  a  panel.  Induced 
velocity  formulas  are  derived  from  the  surface  analogy  of  the  Biot-Savart  Law, 
and  the  edge  vortex  is  avoided.  To  insure  that  the  vorticity  conservation 
theorems  are  satisfied  over  the  interior  of  a  panel,  the  vector  vorticity 
strength  is  expressed  in  terms  of  derivatives  of  an  equivalent  dipole  strength. 

To  the  extent  that  edge  vortices  on  adjacent  panels  do  not  completely  cancel 
the  velocity  field  produced  by  the  sum  of  all  panels  is  nonpotential,  but 
this  error  is  of  order  equal  to  other  errors  of  the  formulation.  To  insure 
this  result,  the  nominal  dipole  strength  must  be  nearly  continuous  across 
panel  edges.  This  is  always  the  case  for  adjacent  panels  in  the  same  lifting 
strip  of  panels  (figure  1).  To  obtain  near  continuity  from  strip  to  strip, 
i.e.,  across  the  N-lines,  it  is  necessary  to  use  the  pi ecewise-1 inear  spanwise 
vorticity  variation  rather  than  the  piecewise  constant  option  (ref.  1).  Along 
physical  edges  of  the  body,  such  as  sheared  off  wing  tips,  an  edge  vortex  is 
"physically"  required.  Either  physical  edges  require  special  formulas  not  incor¬ 
porated  into  the  initial  version  of  the  present  method  or  special  techniques  for 
use  of  the  method  are  needed  to  simulate  edge  effects.  Clearly  further  study 
of  this  situation  is  required  (section  5.3). 

The  principal  analytic  task  of  the  present  work  consisted  of  deriving 
analytic  expressions  for  the  various  integrals  of  section  5.0.  It  is  important 
that  these  induced  velocity  integrals  be  obtained  analytically.  The  alterna¬ 
tive  is  to  perform  numerical  integrations  during  each  production  computer  run 
and  thus  greatly  increase  the  computation  time.  On  the  contrary,  with  the 
integrations  performed  analytically,  the  computing  time  for  the  higher-order 
program  is  about  the  same  as  that  of  the  base  method  for  the  same  number  of 
panels.  (The  former  is,  of  course,  considerably  more  accurate.)  Not  only 
have  general  expressions  for  the  integrals  been  derived,  but  also  simplified 
expressions  to  be  used  to  save  computing  time  when  the  control  point  in 
question  is  far  from  the  panel  (refs.  1  and  2).  The  induced-velocity  formulas 
are  presented  in  section  6.0. 
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Once  the  above  formulas  have  been  obtained,  the  potential  or  velocity 
induced  by  the  source  distribution  on  a  panel  at  an  arbitrary  point  in  space 
(including  any  control  point)  can  be  written  as  the  sum  of  six  terms,  each  of 
which  consists  of  the  product  of  one  of  the  above  integral  expressions  and  a 
certain  coefficient.  The  base-method  integral  has  as  its  coefficient  the 
unknown  value  of  source  density  at  the  control  point  of  the  panel  itself. 

For  three  other  integrals,  each  coefficient  consists  of  the  product  of  this 
source  density  and  one  of  the  three  local  curvatures  of  the  surface.  The 
curvatures  are  evaluated  numerically  by  special  procedures,  which  are  out¬ 
lined  in  section  7.1.  The  curvatures  are  calculated  once  and  for  all  for 
each  panel  before  any  induced- velocity  computation  begins.  Thus,  the  above 
four  integrals  may  be  added  to  obtain  a  single  term  whose  coefficient  is 
the  unknown  value  of  source  density  at  the  control  point  of  the  panel  -  an 
altered  "influence  coefficient"  whose  role  is  exactly  analogous  to  the 
corresponding  coefficient  in  the  base  method. 

The  coefficients  of  the  other  two  integrals  are  the  unknown  derivatives 
of  the  source  density  with  respect  to  distance  in  two  orthogonal  directions 
tangent  to  the  panel.  These  derivatives  must  be  expressed  in  terms  of  the 
unknown  value  of  source  density  at  the  control  point  of  the  panel  and  the 
values  of  source  density  at  the  control  points  of  the  surrounding  panels. 

Once  this  has  been  done,  the  contribution  associated  with  each  surrounding 
value  of  source  density  is  added  to  its  above-described  altered  "influence 
coefficient"  to  obtain  a  final  "influence  coefficient"  that  expresses  the 
total  effect  of  that  value  of  source  density.  This  is  all  very  simple  if 
the  "matrix  of  influence  coefficients"  is  computed  by  rows  or  if  it  can  be 
stored  in  the  high-speed  random-access  core  of  the  computer.  However,  in 
three-dimensional  cases  this  matrix  is  computed  by  columns,  and  it  is  normally 
an  order  of  magnitude  too  large  to  fit  in  core,  so  large  amounts  of  low- speed 
storage  are  required.  It  is  clear  from  Figure  1  that  the  panels  surrounding 
a  particular  panel  consist  not  only  of  panels  on  the  same  strip  but  also 
panels  on  adjacent  strips.  In  the  overall  numbering  scheme  of  the  panels, 
these  latter  may  be  several  tens  of  panels  removed  from  the  panel  in  question. 
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Accordingly,  the  task  of  associating  the  proper  additional  influences  with 
these  panels  is  nontrivial.  A  rather  complicated  logical  scheme  involving  the 
use  of  several  auxiliary  tape  units  has  been  worked  out  and  Implemented.  This 
procedure  is  described  in  section  8.2. 

Similarly  the  induced  velocity  due  to  the  vorticity  distribution  on  a 
panel  is  a  sum  of  terms,  each  of  which  is  the  product  of  a  first  or  second 
derivative  of  the  equivalent  dipole  distribution  (value  and  first  derivatives 
of  vorticity)  and  some  integral  expression.  Some  terms  also  contain  the  local 
geometric  curvatures.  The  above  elaborate  procedure  that  accounts  for  the 
source  derivative  effects  is  required  because  the  source  is  the  basic  singu¬ 
larity.  Since  the  vorticity  is  an  auxiliary  singularity,  simpler  procedures 
may  be  used  to  account  for  its  derivatives.  The  one  adopted  in  the  initial 
version  of  the  present  method,  which  is  described  in  section  9.0,  is  the 
simplest  of  all.  It  could  be  made  more  elaborate  and  presumably  more  accurate 
without  greatly  complicating  the  procedure. 
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5.0  CONSISTENT  EXPANSIONS  FOR  THE  POTENTIAL  AND  VELOCITY  INDUCED 
BY  A  CURVED  PANEL  AT  A  POINT  IN  SPACE 


5.1  The  Potential  of  a  Source  Distribution 

A  "true"  panel  on  a  body  surface  is  the  curved  four- sided  region  of  the 
surface  whose  corners  are  input  points  lying  exactly  on  the  surface  (figure  3). 
The  boundary  curves  of  the  true  panel  that  connect  the  corner  points  will  be 
defined  shortly.  Consider  a  plane  tangent  to  the  surface  at  some  central 
point  of  the  true  panel  as  yet  unspecified.  A  panel  coordinate  system  is  con¬ 
structed  whose  origin  is  the  tangency  point  and  whose  z  or  c  axis  is  normal 
to  the  tangent  plane.  The  x  or  c  and  y  or  n  axes  lie  in  this  plane,  and 
their  precise  orientation  is  unimportant  for  purposes  of  this  section.  The  cor¬ 
ner  points  of  the  true  panel  are  projected  into  the  tangent  plane  along  the  nor¬ 
mal  direction.  By  joining  adjacent  projected  corner  points  with  straight  lines, 
a  flat  panel  is  produced,  which  is  assumed  to  be  the  projection  of  the  true  panel 
in  the  tangent  plane.  This  construction  now  defines  the  boundary  curves  of  the 
true  panel.  They  are  the  curves  joining  the  true  corner  points  that  have 
straight-line  projections  in  the  tangent  plane. 


It  is  a  fundamental  assumption  of  the  present  method  that  the  dimensions 
of  the  panel  are  small  in  certain  senses.  Certainly  variations  over  the  panel 
of  the  normal  direction  are  assumed  small.  Moreover,  if  a  physical  quantity 
is  expanded  in  a  Maclaurin  series,  successive  terms  become  small,  i.e.,  if 


00  00 

(s  3^  +  n  a^r) o,o 

n=0  '  n=0 


then 


IQJ  «  IQJ 


if 


m  <  n 


0) 

(2) 


Also,  the  vertical  distance  c  of  a  point  on  the  true  panel  is  of  the  order 
of  the  square  of  the  horizontal  distance,  i.e., 

t  =  OU2  +  n2)  «4>k*  +  n2  (3) 
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The  potential  at  a  point  (x,y,z)  in  space  due  to  a  source  distribution 
on  the  true  panel  is 


♦  ■  ffvi%  <4> 

S 

where 

r2  =  (x  -O2  +  (y  -n)2  +  (z  -  c)2  (5) 

and  w'uere  S  is  the  surface  area  of  the  true  panel.  It  is  desired  to  express 
<t>  in  terms  of  a  series  of  integrals  over  the  projected  flat  panel.  Thus  it 
may  be  stated  that  the  potential  of  the  true  panel  is  "expanded  about"  that 
of  the  flat  panel.  To  illustrate  the  process  more  completely,  a  three-term 
expansion  is  derived,  although  only  the  first  two  terms  are  ultimately 
retained  in  the  present  method. 

In  panel  coordinates  the  equation  c  =  f(£,n)  of  the  surface  of  the 
true  panel  may  be  expanded  in  a  Maclaurin  series  in  the  form 

C  =  [P£2  +  2QCn  +  Rn2]  +  [T3q£3  +  T2^£2n  +  T125n^  +  T03n3^  +  (6) 

=  C2  +  c3  +  ... 

There  are  no  constant  or  linear  terms  in  (6)  because  the  origin  is  at  the 
tangency  point.  All  coefficients  in  (6)  are  constants  proportional  to 
derivatives  of  c  at  the  origin.  The  coefficients  P,  Q,  and  R,  which  are 
the  only  ones  actually  used  in  the  present  method,  are  the  second  derivatives. 
They  are  referred  to  below  as  the  surface  curvatures  to  which  they  are  closely 
related.  The  quantity  ?2  1S  second  order  in  £  and  n  and  thus  in  panel 
dimension  t,  and  c3  is  third  order. 

The  equation  of  the  true  panel  may  be  written 

FU.n.c)  =  c  -  S2(c»n)  “  53(Cin)  ~  ...  *  0  (7) 

Then,  taking  the  gradient  gives 

grad  F  =  t  -  U2?  +  +  ...)i  -  (c2n  +  C3n  +  «..)j  (8) 

where  subscripts  £  and  n  denote  partial  derivatives  and  T,  j.  It  are  the 
unit  vectors  of  the  panel  coordinate  system.  The  vector  grad  F  is  normal 
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to  the  panel  at  any  point.  The  unit  normal  at  any  point  is 


n  =  grad  F 
"[grad  F| 

so  the  c-component  of  the  unit  normal  is 

1  _  1 


"c  = 


I  grad  F| 


1  +  +  C3c  +  *•*)  +  (t2n  +  C3n  + 


(9) 


00) 


This  can  be  expanded  in  the  form 


nC  =  1  2  ^c2c  +  c3? 


+  (c2n  +  c3n 


*  8  +  C3C  +  •••>  +  (;2„  *  s3n  *  +  •••  (11) 

From  (6)  it  is  clear  that  ^  and  ^  are  first  order,  and  $3  are 

second  order,  etc.,  so  the  leading  terms  from  the  first  square  bracket  of  (11) 
are  second  order,  and  the  leading  terms  in  the  second  square  bracket  are 
fourth  order.  Thus, 


n 


c 


(12) 


is  a  valid  three-term  expansion  with  second  (linear)  term  zero.  More  pre¬ 
cisely, 

n?  =  1  -  \  [(2PC  +  2Qn)2  +  (2QC  +  2Rn)2] 

=  1-2  [(P2  +  qV  +  2(PQ  +  QR)Cn  +  (Q2  +  R2)n2] 

n?  =  1  -  2ip2 

where  ip2  second  order.  The  elementary  surface  area  dS  on  the  true 
panel  is  related  to  the  elementary  area  dA  =  dedn  in  the  tangent  plane  by 


n  dS  =  dA  (14) 

dS  =  (1  +  2<j>2)d£dn  (15) 

The  source  density  is  strictly  a  function  of  surface  distances  along  the  panel. 
However,  it  can  be  shown  that  there  is  no  difference  between  surface  distances 
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and  e  and  n  through  second  order, 
of  5  and  n  in  the  form 

o  =  ao  +  (ax£  +  ayn)  +  (o^2 


Thus  it  suffices  to  define 
+  2axyen  +  oyyn2)  +  ... 


a 


in  terms 


=  cr  +  a,  +  a0 
0  1  c. 

where  on  is  n-th  order  in  £  and  n  and  thus  in  panel  dimension 
coefficients  in  (16)  are  constants  proportional  to  derivatives  of  0 
origin. 


(16) 

t.  The 
at  the 


Thus,  to  three  terms 


where 


0dS  *  (aQ  +  a1  +  o2)(l  +  2^2)d£dn 
=  (°0  +  +  apdCdn 


(17) 


a2  =  a2  +  2^2  (18) 

All  of  the  above  expansions  are  independent  of  the  location  of  the  point 
(x,y,z). 


It  remains  to  expand  1/r,  which  of  course,  does  depend  on  x,y,z.  It 
is  necessary  to  differentiate  three  ranges  of  value  of  r: 


(a) 

r  »  t 

for  all  c.n 

(b) 

-5 

II 

O 

for  all  £,n 

(19) 

(c) 

r  «  t 

r  =  0(t) 

for  some  5,n 

for  other  $,n 

In  words,  the  ranges  amount  to  the  situations  where  the  distance  of  the  point 

(x.y.z)  from  the  origin  of  panel  coordinates  is,  respectively  large,  of  the 

order  of,  and  small  compared  to  the  dimensions  of  the  panel.  It  turns  out 

that  the  range  (b)  where  r  is  of  the  order  of  the  panel  dimensions  is  the 

essential  one.  In  the  far-field,  range  (a),  1/r  is  expanded  in  negative 
2  2  2  h 

powers  of  rQ  =  (x  +  y  +  z  )  .  The  resulting  expansion  differs  in  no 
important  way  from  the  usual  far-field  multipole  expansion  and  thus  is 
relatively  easy  to  derive.  An  order-of-magnitude  comparison  of  the  expansions 
of  ranges  (a)  and  (b)  shows  that  while  certain  terms  become  small  faster 
than  others  as  the  distance  r  increases,  the  expansion  derived  for  range  (b) 
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is  valid  for  range  (a),  although  it  is  conservative  in  that  some  retained 
terms  could  be  eliminated.  For  points  close  to  the  panel,  range  (c),  it  is 
necessary  to  assume  that  they  lie  on  a  line  through  the  origin  having  a  finite 
slope  with  respect  to  the  tangent  plane.  Under  this  condition  the  order-of- 
magnitude  analysis  of  range  (b)  remains  valid.  This  is  just  what  the  physics 
of  a  panel  method  requires  in  any  event.  Eventually,  the  control  point  of 
the  panel  is  identified  with  the  origin  of  panel  coordinates  and  the  above 
condition  states  that  if  a  point  approaches  the  surface,  it  does  so  at  a 
control  point.  Since  it  is  only  at  the  control  points  that  the  normal  velocity 
boundary  condition  is  applied,  approaching  any  other  surface  point  would  give 
physically  meaningless  results.  Thus  the  derivation  below  concentrates  on 
range  (b),  which  may  be  thought  of  as  the  effect  of  a  panel  on  control  points 
of  nearby  panels. 

The  distance  r  can  be  written 

r2  =  (x  -  O2  +  (y  -  n)2  +  z2  -  2zs  +  x,2 

=  r2  -  2zc  +  c2  (20) 

where  rf  is  the  distance  between  (x,y,z)  and  the  point  (c.n.O)  on  the 
flat  panel  (Figure  3).  Thus 


Note  that  (z/rf)  is  of  order  unity,  and,  since  r^  is  0(t),  the  quantity 
c/rf  is  also  0(t),  i.e.,  0(c).  If  use  is  made  of  (6),  the  above  becomes 


from  which  the  desired  three-term  expansion  is  obtained  in  the  form 
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or,  for  abbreviation. 


—  =  - —  (1  +  C,  +  Co) 
r  rf  v  1  2 ' 

where  cn  is  of  order  n  in  £  or  n  and  thus  in  t. 
by  (17)  to  obtain  the  expansion 

r  dS  =  r^"  ^ao  +  °1  +  a2^1  +  C1  +  c2^dA 

=  7^  Ca0  +  CT1  +  °2 

+  cl°o  +  clal  +  cl°2 


(23) 

Now  this  is  multiplied 


(24) 


+  C2a0  +  C2a]  +  C^] 

The  square  bracket  in  (24)  must  be  reduced  to  a  three-term  expansion  using  the 
facts  that 


a  »  a i  >>  Oo 
o  I  c. 

and 


(25) 


1  »  c-j  »  c2  (26) 

The  leading  (lowest  order)  term  is  clearly  aQ  because  all  terms  are  small 
compared  to  it.  Possible  members  of  the  second  term  are  and  c-jo0 
because  all  remaining  terms  are  small  compared  to  one  or  the  other.  Thus, 
both  of  these  are  retained  in  the  second  term,  because  neither  can  be  guaran¬ 
teed  to  be  small  compared  to  the  other  in  all  cases.  The  question  then  arises 
as  to  which  of  the  remaining  six  terms  of  the  square  bracket  of  (24)  need  to 
be  retained  in  the  third  term  of  the  expansion.  Obviously,  if  any  of  these 
six  is  small  compared  to  any  other,  it  may  be  discarded.  This  eliminates  all 
but  o2,  Ci o i ,  and  c2oQ.  Thus,  the  three-term  expansion  of  the  integrand 
of  (4)  has  the  form 

r  dS  =  7^  Cao  +  (cl°o  +  Gl)  +  (c2ao  +  cl°l  +  a2)]  (27) 

when  the  abbreviations  of  (16)  and  (24)  are  replaced  by  their  actual  expres¬ 
sions,  the  three-term  expansion  of  (4)  is 


+  [a0  ff( (Z  +  ( V  +  ayn)dA  (28> 

L  A  '  f  f  rf 7  A  rf 

+ff  T7  \axx^  +  2cxy^  +  V"2  +  2(p2  +  q2)??  +  4(PQ  +  QR)Cn 

A 

+  2(Q2  +  R2)n2  dA 

where  the  integrals  are  over  the  projected  flat  panel.  Defining  the  integrals 

1  =  ff  -^TT  dA  (29) 

mnp  JJ  p 

A  rf 

the  three-term  expansion  of  the  potential  can  be  written 
*  =  </°>a  +  L(c)a  +  *(lx)a  +  «f,(ly>o 

0  1  0  x  y  (30) 

+  <t>^20^a  +  4>^2x^0  +  </2y)a  +  4>^2xx  )<j'  +  2<t>(2xy)a'  +  ^2yyK' 

where  i  0  x  y  **  xy  yyj 


o'  =  a  +  2(P2  +  021 
XX  xx  vr  *  > 


“Jy  =  °xy  *  2W  *  Q«) 
°'yy  =  °yy  +  2«2  +  r2> 


The  individual  potentials  in  (30)  are 


4’(  *  "  z^T30I303  +  T21 X213  +  T12I123  +  T03!033^ 

+  I  Z'  tP2r405  +  4P^315  +  <2PR  +  4^)I225  +  4(>RI135  + 

-  1  [p2l403  +  4PQI313  +  (2PR  +  4(52)I223  +  4QRI133  +  r2i043]  <35> 

<fi(  )  =  z[PI3Q3  +  2QI213  +  RI123] 

(P  )  (36) 

*{' ^  =  Z[PI213  +  ^123  +  RW 


$(2xx)  =  I, 

t 

/2*y> .  i. 
»(2yy)  .  - 


(37) 


The  first  term  of  (30),  ,  corresponds  to  a  flat  panel  with  a 

constant  source  density.  This  is,  of  course,  the  term  used  in  the  base  method. 

The  second  term  of  (30)  contains  the  second  derivatives  P,  Q,  and  R,  but 
no  higher  derivatives,  of  the  surface  shape  and  contains  the  first  derivatives 

of  the  source  density,  but  no  higher  derivatives.  Thus  the  second  term  of  (30) 

corresponds  to  a  paraboloidal  panel  shape  with  a  linearly  varying  source  den¬ 
sity.  The  third  term  of  (30)  contains  all  the  preceeding  quantities  and  also 
the  third  derivatives  of  the  surface  shape  and  the  second  derivatives  of  the 
source  density:  a  cubic  panel  with  a  quadratic  source  density. 


The  equations  above  illustrate  the  fact  that  succeeding  terms  in  the 
expansion  for  the  potential  of  a  panel  increase  rapidly  in  complexity.  A 
single  first  term  is  followed  by  a  second  term  containing  five  individual 
parts,  each  with  its  own  integral  of  the  form  (29).  The  third  term  contains 
23  individual  parts  which  together  involve  17  different  integrals  of  the  form 
(29).  The  great  increase  in  complexity  associated  with  retaining  the  third 
term  of  (30)  appears  to  be  unjustified  at  this  time.  Accordingly,  the  higher- 
order  method  accounts  for  the  source  density  effect  by  considering  the  first 
two  terms  of  (30).  This  is  the  approach  that  has  previously  been  followed  in 
the  two-dimensional  and  axisymmetric  higher-order  methods,  references  5  and  6. 
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5.2  The  Potential  of  a  Dipole  Distribution 

This  parallels  almost  exactly  the  development  of  section  5.1  for  the 
effect  of  a  source  distribution  on  a  panel,  as  sketched  in  figure  3.  If  the 
dipole  strength  is  p(E,n)  the  potential  of  the  true  panel  is 


♦d  ■  jO1]"  •  9rad£nC  <?>|  ■  fft*c 


where  as  before  n  is  the  unit  normal  vector  to  the  true  panel.  From  equa¬ 
tions  (7)  through  (13)  it  is  known  that  (13)  is  a  valid  three  term  expansion 
of  both  the  c-component  of  n  and  of  |gradF|~\  Expansions  of  the  other 
two  components  of  the  normal  vector  are 

"5  m  -(c2t  +  '3c)ne  m 

\  -  -<e2r,  +  Wt 

Only  two  terms  are  required  because  these  two  are  respectively  of  first  and 
second  order  in  e  and  n  while  n^  has  a  zero  order  term.  Thus 

n  •  grad(^)  =  [-(s2  +  “  O  -  (^n  +  c3n^y  “  +  z  ~  ^  (40) 


is  a  valid  three-term  expansion,  and  it  is  clear  retention  of  third  terms  in 
(39)  would  produce  fourth  terms  in  (40).  When  use  is  made  of  (14)  the  inte¬ 
grand  of  (38)  becomes 


d<f>d  =  {z  -  [c2  +  ?2e(x  -  e)  +  ?2r)(y  -  n)] 


-  [c,  +  ?-jr(x  -  C)  +  Co  (y  -  n) ] } 


p0  +  vi  +  w2 


where 

pQ  =  const 

“i  =  NE  +  V  < 

2  2 
p_  =  u  e  +  2u  En  +  p  n 
*2  Mxx*  Hxy*  1  Hyy 1 

and  where  the  terms  in  the  curly  bracket  have  been  arranged  in  increasing 
order. 
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From  (23) 


■T  =  T’  (1  +  C1  +  c2)3 

r  r^ 

•  ^  [1  +  3c]  +  3(c^  +  c2)] 
rf 


(43) 


By  reasoning  similar  to  that  following  equation  (26)  a  valid  three-term 
expansion  is 


Finally  a  valid  three-term  expansion  of  (41)  is 

d*d  =  {Zwo  +  {z(3clvo  +  “  [?2  +  c2?(x  “  +  c2n(y  “  n)]uo} 

+  {z[3(c^  +  c2)v0  +  30^-j  +  v2)  “  I?2  +  ?2^x  “  +  C 2r) (y  ~  n)](3c1u0  + 

~h  +  WX  "  C)  +  c3n(y  "  n)  V  dA  (45) 


The  problem  with  using  dipoles  is  that  the  three-term  expansion  (45)  of 
the  dipole  potential  (38)  is  equivalent  in  accuracy  to  a  two-term  expansion 
of  the  source  potential  (4)  obtained  by  truncating  (30).  This  can  be  verified 
in  several  ways.  Mathematically  it  is  straightforward  to  verify  that  when 
x,y,z  and  r^  are  of  order  5  and  n  and  thus  of  the  order  of  element 
dimension  t,  each  term  of  (45)  is  one  order  lower  than  the  corresponding 
term  of  (30).  For  example  the  leading  terms  are 


Vf  = 

Vf  *  -o011'2' 


(46) 


Some  physical  insight  can  be  obtained  by  performing  the  integration  for  the 
first  term  of  (45)  which  corresponds  to  a  constant-strength  dipole  distribu¬ 
tion.  The  result  is  the  potential  of  a  concentrated  vortex  filament  lying 
around  the  edges  of  the  panel.  Clearly  no  such  singularity  can  exist  on  the 
body  surface.  It  is  present  only  because  an  isolated  panel  has  been  con¬ 
sidered,  and  it  must  cancel  (at  least  to  some  order)  with  similar  singular¬ 
ities  on  the  surrounding  panels.  Thus,  when  only  noncancelling  terms  are 
considered,  a  valid  two- term  expansion  of  (38)  consists  of  the  second  and 
third  terms  of  (45),  which,  as  mentioned  above,  correspond  in  their  orders 
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to  the  first  two  terms  of  (30).  Unfortunately  they  are  much  more  complicated. 
Integration  of  the  third  term  of  (45)  leads  to  an  expression  even  more  compli¬ 
cated  than  the  third  term  of  (30).  Another  difficulty  is  that  the  third 
derivatives  of  the  surface  shape  enter  into  the  third  term  of  (45)  through 
and  its  derivatives.  Thus,  although  only  quadratic  dipole  terms,  which 
correspond  to  linear  vorticity  terms,  are  included  in  (45),  it  would  appear 
that  a  cubic  panel  is  required.  However,  this  turns  out  to  be  only  apparently 
true. 


Most  of  the  difficulty  with  the  expansion  (45)  arises  from  the  fact  that 
it  expresses  the  effect  not  only  of  a  distributed  vorticity  distribution  over 
the  panel,  but  also  of  a  concentrated  vortex  filament  around  the  panel  edges. 
As  mentioned  above,  the  first  term  represents  a  constant-strength  vortex 
filament  and  no  distributed  vorticity  at  all.  The  second  and  third  terms 
contain  mixtures  of  effects  due  to  distributed  vorticity  and  those  due  to 
variable-strength  vortex  filaments  around  the  panel  edges.  In  particular, 
all  the  third  derivative  terms  enter  only  by  way  of  these  "edge  vortex"  terms. 
Clearly  all  edge  vortex  terms  must  cancel  with  those  of  adjacent  panels.  The 
result  must  be  a  much-simplified  expansion  and  one  that  does  not  contain  any 
third  derivative  terms.  Unfortunately,  it  does  not  appear  possible  to  accomp¬ 
lish  this  cancellation  analytically.  Certainly  it  is  necessary  to  use  the 
integrated  form  of  (45),  because  the  cancellation  is  not  an  identity  in  the 
integrand,  but  only  in  the  final  result.  Even  if  this  is  done,  the  effects 
of  the  distributed  vorticity  and  those  of  the  edge  vortex  are  not  neatly 
separated  but  are  combined  into  some  of  the  analytic  expressions.  Clearly 
all  terms  containing  yQ  should  cancel,  but  it  is  not  sufficient  simply  to 
drop  these  from  (45).  Unless  all  cancelling  terms  can  be  identified  and 
eliminated,  the  deletion  of  some  of  them  leads  to  a  disruption  of  the  cancel¬ 
lation  process. 

5.3  The  Velocity  Due  to  a  Vorticity  Distribution 

The  situation  described  in  section  5.2  is  present  also  in  two  dimensions. 
Indeed,  the  considerations  of  the  previous  section  regarding  cancellation, 
ecc.,  can  be  illustrated  most  easily  with  two-dimensional  examples.  Dif¬ 
ficulties  are  avoided  in  two  dimensions  by  considering  a  vorticity  distribu¬ 
tion  on  the  panel  rather  than  a  dipole  distribution.  In  such  a  case  all  the 
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consistency  arguments  of  sections. 1  can  be  invoked,  because  the  velocity  due 
to  vorticity  is  simply  that  due  to  the  same  distribution  of  source  strength 
rotated  through  a  right  angle.  This  approach  simply  ignores  edge  vortex  terms 
on  the  presumption  that  they  cancel  with  those  of  adjacent  panels. 

The  present  method  applies  this  technique  to  the  three-dimensional  case, 
i.e.,  it  uses  a  linearly-varying  vorticity  distribution  over  the  panel  and 
ignores  edge  vortex  effects.  (Reference  1  gives  specific  proof  of  the  rela¬ 
tion  between  dipole  and  vorticity  distributions  and  the  edge  vortex.)  There 
are,  of  course,  some  differences  between  the  two-dimensional  and  the  three- 
dimensional  cases.  The  most  obvious  theoretical  danger  arises  from  the  fact  that 
the  velocity  field  due  to  a  vorticity  distribution  is  in  general  nonpotential 
with  edge  vorticies  neglected.  To  the  extent  that  edge  vortex  cancellation  with 
adjacent  panels  is  less  than  exact,  the  final  velocity  field  about  the  body  is 
nonpotential.  However,  this  does  not  appear  to  be  serious.  Any  panel  method  is 
an  approximation  to  the  exact  potential  flow  with  calculational  errors  of  some 
higher  order  that  approach  zero  with  increasing  panel  number.  There  is  no  obvious 
reason  why  a  method  whose  calculation  error  is  potential  should  be  preferred  to 
one  whose  error  is  not.  There  are,  however,  two  features  of  the  use  of  vorticity 
worth  mentioning.  Surface  vorticity  is  an  auxiliary  singularity  in  the  sense  of 
section  3.0.  While  its  variation  over  a  single  panel  is  linear,  its  global 
variation  over  the  body  as  a  whole  is  legislated  by  what  may  be  termed 
"assembly  formulas"  (section  9.0)  that  contain  just  enough  free  parameters  to 
permit  satisfaction  of  the  Kutta  condition  on  the  trailing-edge  segments 
(figure  1).  It  is  important  that  the  assembly  formulas  insure  approximate 
continuity  of  the  dipole  strength  that  is  equivalent  to  the  surface  vorticity 
so  that  edge  vortex  effects,  in  fact,  approximately  cancel.  In  particular,  the 
step- function  spanwise  variation  of  bound  vorticity  (ref.  1),  which  is 
ordinarily  used  with  the  base  method,  is  not  appropriate  for  use  with  the 
higher-order  surface  vorticity.  Finally  it  is  clear  that  at  any  actual  edges 
of  the  body  where  a  concentrated  vortex  filament  is  truly  required,  special 
edge- vortex  formulas  must  be  incorporated  into  the  method.  The  most  common 
such  situation  is  a  wing  tip  with  a  nonzero  chord,  if  the  equivalent  dipole 
strength  is  not  forced  to  go  to  zero  there.  In  practical  cases  only  a  small 
percentage  of  the  panels  would  require  these  additional  formulas.  However, 
at  this  time  edge-vortex  formulas  have  not  been  included,  so  that  the  method 


25 


is  not  applicable  to  any  body  requiring  edge  vortices.  Inclusion  of  edge- 
vortex  formulas  is  seemingly  a  desirable  modification  of  the  method.  In 
the  meantime,  however,  cases  requiring  edge  vortices  must  be  avoided.  For 
example,  a  wing-tip  is  handled  either  by  closing  the  tip  to  a  point  or  by 
having  the  assembly  formulas  force  the  equivalent  dipole  strength  to  zero  at 
the  tip.  This  last  can  be  accomplished  in  a  routine  manner  using  the  piecewise- 
linear  spanwise  vorticity  variation  (section  9.4).  Either  or  both  of  these 
procedures  should  give  satisfactory  results  for  a  normal  wing  tip.  However, 
a  partial-span  deflected  flap  may  be  more  difficult  and  requires  further 
study. 

The  effect  of  a  vorticity  distribution  on  a  panel  cannot  be  expressed  in 
terms  of  a  potential,  as  mentioned  above.  The  velocity  induced  by  the  panel 
of  figure  3  at  a  point  (x,y,z)  is 


where  <o  is  the  vector  vorticity  strength  and 

r  =  (x  —  5)T  +  (y  -  n)j  +(z  —  ?)k  (48) 

The  distance  r  has  its  usual  meaning  (5),  which  is  also  equal  to  |r|,  and 
the  integral  is  over  the  true  panel.  To  insure  that  the  vorticity  satisfies 
the  usual  vorticity  conservation  theorems  over  the  panel,  it  is  convenient  to 
express  u  in  terms  of  an  equivalent  dipole- distribution  y.  As  shown  in 
reference  1,  the  relation  is 

u>  =  -n  x  grad  y  (49) 

where  n  is  the  unit  normal  vector  (9),  whose  components  n^»  are 

given  by  (13)  and  (39).  To  be  compatible  with  the  source  density  and  panel- 
geometry  expansions,  y  is  taken  to  be  a  quadratic  function  of  panel  coord¬ 
inates  ?  and  n  in 

v  =  vo  +  px5  +  V  +  pxx*2  +  lJxy£n  +  pyyn2  (5°) 

so  that  the  components  of  m  vary  linearly  over  the  panel.  Furthermore, 
since  a  two- term  expansion  of  the  source  potential  (4)  is  all  that  appears 
feasible,  only  a  two-term  expansion  of  (47)  is  required.  Evidently 
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=  p  +  2(p  5  +  p  n) 
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are  two-term  expansions  of  the  components  of  gradp.  Then  (49)  gives 
u  =  (n  iE)T  +  (n  in  -  n  in\£ 

vnt  3n*  v  C  85' J  inn  85  5  8n' 


(52) 


Two-term  expansions  of  the  i  and  j  components  of  (52)  contain  zeroth  and 
first-order  terms.  Since  the  leading  term  of  the  k  component  is  linear  only 
that  one  term  is  required.  The  two-term  expansion  of  (52)  is 

21  ’  lvy  *  2("xy{  +  V>))T  ~  <\  +  2("xxe2  +  V2)1^  *  f-'2„»x  +  E25V* 


=  co_i  +  a)  j  +  a)  k 
S  n  c 

Also  needed  is 


(53) 
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For  the  present  purpose  the  expansion  (43)  for  r~^  may  be  truncated  after 
two  terms,  i.e. , 


(55) 


Moreover  to  the  order  being  considered  n?*l,  and  thus  dS  =  dA  (15).  If 
now  the  components  of  (47)  are  written 
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the  integrands  are 
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dV 
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where 
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*z  -  (x  -  c)px  +  (y  -  n)uy 
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(57) 


(58) 


Jz  =  2^pxx5  +  uxyn^x  “  ^  +  ^pxy5  +  wyyn^y  “ 

In  each  of  the  three  components  of  (57),  the  first  term  is  the  leading  onewitha 
o 

magnitude  0(t“  )  neglecting  the  factor  dA.  These  leading  terms  correspond 
to  a  constant  vorticity  on  a  flat  panel.  The  second  and  third  terms  in  each 
component  of  (57)  have  the  same  magnitude  0(t"^)  and  together  comprise  the 
higher-order  term  of  the  expansion.  The  third  term  is  a  pure  curvature  effect, 
while  the  second  term  has  both  curvature  and  vorticity  derivative  effects. 

Regardless  of  their  physical  nature,  the  second  terms  of  (57)  have  the 
mathematical  form  of  potentials  of  linear  and  quadratic  dipoles  on  a  flat 
panel,  and  thus  they  can  be  evaluated  by  the  formulas  of  reference  1.  The 
same  is  true  of  the  first  terms.  Only  the  third  terms  of  (57)  require  new 
analytic  formulas. 


The  above  considerations  make  it  convenient  to  consider  the  vortex 
velocity  as  a  sum  of  the  three  terms  in  the  form 
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wZ 

UZ 
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(59) 


•where  the  superscripts  0,  1,  2  denote,  respectively,  integrals  of  the  first, 
second,  and  third  terms  of  (57). 
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6.0  INTEGRATED  INDUCED  VELOCITY  FORMULAS 
FOR  A  PANEL 


6.1  General  Remarks 

This  section  presents  detailed  formulas  for  the  velocity  induced  by  an 
individual  panel  at  a  point  in  space.  These  formulas  are  obtained  from  the 
expressions  developed  in  section  5.0  by  performing  the  indicated  Integrations 
over  the  flat  projected  panel.  Different  procedures  are  called  for  depending 
on  the  distance  of  the  point  in  question  from  the  panel.  For  nearby  points 
the  expressions  of  section  5.0  are  integrated  exactly.  This  procedure  is 
rather  lengthy  and  details  are  omitted.  Only  the  final  formulas  for  this 
"Near  Field"  are  presented,  and  these  are  the  key  to  the  present  method.  As 
in  reference  1,  it  is  assumed  that  the  point  in  question  has  been  transformed 
into  the  panel  coordinate  system  and  all  Near-Field  formulas  are  given  in 
terms  of  this  coordinate  system.  For  points  further  from  the  panel,  the 
integrals  of  section  5.0  are  evaluated  by  a  classic  multipole  expansion 
(refs.  1  and  2).  The  orders  of  the  expansions  are  selected  to  be  at 
least  as  high  as  the  terms  in  question.  This  is  a  relatively  simple  procedure 
analytically,  and  the  resulting  "Intermediate  Field"  formulas  require  much 
less  computing  time  than  the  Near-Field  formulas.  This  computation  also  is 
carried  out  in  panel  coordinates.  For  points  even  further  away,  a  "Far-Field" 
approximation  is  used.  This  is  obtained  simply  by  retaining  only  the  first 
terms  in  the  multipole  expansions.  However,  for  the  source  effects  <>  d  some 
of  the  vorticity  effects,  the  Far-Field  formulas  have  been  put  in  vector  form, 
and  thus  they  can  be  evaluated  directly  in  the  so-called  reference  coordinate 
system  in  which  the  body  is  input.  This  eliminates  the  need  for  transforma¬ 
tions  and  further  reduces  computing  time.  Certain  terms  in  the  vorticity 
effects  have  such  complicated  Far-Field  expressions  that  it  did  not  seem 
worthwhile  to  put  them  in  vector  form.  Accordingly,  they  are  evaluated  in 
panel  coordinates  as  is  done  in  the  Near-  and  Intermediate-Fields.  Thus,  in 
section  6.0  all  individual  coordinates  and  vector  components  are  in  panel 
coordinates.  Expressions  in  vector  form  may  be  evaluated  in  the  reference 
coordinate  system. 

The  entire  philosophy  of  the  induced  velocity  computation  and  all  of  the 
associated  bookkeeping  details  are  identical  to  those  of  reference  1,  to  which 
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the  reader  is  referred.  The  procedure  for  the  multi  pole  expansion  has  already 
been  mentioned.  Other  matters  are  the  ranges  of  distance  adopted  for  Near-, 
Intermediate-  and  Far-Fields,  transformations  between  panel  and  reference 
coordinate  systems,  and  geometric  quantities  defining  the  projected  flat 
panels.  In  addition,  many  of  the  induced- velocity  expressions  developed  in 
reference  1  enter  the  formulas  of  the  present  method.  To  report  these  would 
add  greatly  to  the  length  of  the  present  report,  so  again  reference  1  is  relied 
on.  The  panel  is  assumed  to  have  two  parallel  sides,  as  illustrated  in  figure  20 
of  reference  1 . 


6.2  Source  Velocity  Formulas 

The  source  potential  is  given  in  equations  (30)  and  (33)  of  section  5.1. 
The  velocity  induced  by  the  panel  is  obtained  by  taking  the  negative  gradient 
to  obtain 

+  IP^  +  2Q^Q)  +  R^R)]  +  ?0x)oy  +  (60) 

o  x  y 

where  each  individual  velocity  is  the  negative  gradient  of  the  corresponding 
potential. 


6.2.1  Near  Field 

The  velocity  components  are: 


V<°>  =  V  (source) 

A  A 

V<°>  =  Vy  (source) 
V<°>  =  V2  (source)  , 


eq.  (7.7.7),  reference  1 
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eqs.  (7.7.14)  and  (7.7.17) 
reference  1  '  (62) 


* 


30 


»; 

V1 

) 

v« 


f(p)  =  _! 

3^20 

z  ~~  +  2x 

X  1 

ax 

r(P) 
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where  all  derivatives  except  JgQ  are  given  in  equations  (7.7.9)  through 
(7.7.11)  of  reference  1,  and  where 


J20  =  (source)  -  H, 


'02 


(64) 
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-Vx (source)  -  3x 
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aj 
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3H, 


ay  =  -Vy (source) 
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3J20  3Hn? 

=  -V2 (source)  -  -  02 


az 


(65) 


Hq2  and  its  derivatives  are  given  in  equations  (7.7.15)  and  (7.7.16)  of 
reference  1 ,  and 


♦ (source 


)  =  Jy-,^.02)  .  (X  ~  {2>  -  "32<y  ~  "2>  L 


(23) 


VT7Z 
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_  (y  _n  )L(34)  -  (_X  -  C4}  m41(y  n4)  l(41)  _ 


(66) 


Vl  +  ml 


zVz (source) 


The  undefined  quantities  in  (66)  are  given  in  section  7.2  and  equation  (7.7.3) 
of  reference  1. 


V<lx)  =  xVx  (source)  -  J2Q 
Vy1X^  =  xVy  (source)  - 


(67) 


,Ox) 


xVz  (source)  -  zVx  (source) 


31 


v^>  =  yVx  (source)  - 
vy1y^  »  yVy  (source)  -  J02 

( 1  w  \ 

J  '  =  yV2  (source)  -  zVy  (source) 


J02  =  H02  “  zVz  (source) 


6.2.2  Intermediate  Field 

Unlike  the  near-field,  where  many  of  the  definitions  of  reference  1  can 
be  used  without  change,  the  intermediate-field  formulas  must  be  modified  by 
the  inclusion  of  terms  involving  first  moments  of  the  panel  area.  These  are 
no  longer  zero  because  the  panel  centroid  is  no  longer  used  as  origin. 

Thus,  for  example,  equation  (7.6.7)  of  reference  1  is  modified  to 


7[->00ux 


xx+I01uxy]  +  (rj  [I20uxxx  +  2I11uxxy 

+  I02uxyy]l 


Vy0)  =  7  j’Wy  +  (r  )[I10uxy  +  r01uyy]  +  (r  )  [I20uxxy  +  2I11U> 


=  7  |_I00uz  +  (^o)[I10u 


V(R)  -  _Z02 

X  »X  ' 

V(R)  _  a*02  , 

y  ay 

V(R)  -  _!!o2 

Z  3Z  * 


+  I_  -  u  1  + 
xz  01  yzJ 


,(Q)  «  _!*n 

x  ax 

/ ( q)  _  _!!n 


(q)  _  3*n 


+  I02uyyy1 


[I20uxxz  +  2IllLxyz 

+  I02‘VyzIi 


eqs.  (7.6.11)  and 
(7.6.12),  reference  1 


It  should  be  noted  that  early  versions  of  reference  1  were  missing  a  factor 
of  1/2  in  the  2nd  order  terms  of  these  equations.  It  should  be  further  noted 
that  when  the  derivatives  of  the  flat-panel  dipole  potentials,  <t>0Q»  $qi  and 


<j>10  are  being  evaluated  in  the  intermediate  field  from  equations  (7.6.8), 
(7.6.9)  and  (7.6.10)  of  reference  1,  not  only  should  the  second-order  terms 
be  multiplied  by  a  1/2  (if  such  a  factor  is  not  already  present),  but  also 
all  "crossed  out"  terms  should  be  included.  This  last  is  due  to  the  fact  that 
the  panel  centroid  is  no  longer  used  as  origin  of  panel  coordinates  and  thus 
the  first  moments  1^  and  are  not  zero. 

(P)  _  a*20 

X  3X 

(72) 

y  3y  v  ' 

(P)  =  a*20 

Z  3Z 

where  and  its  derivatives  are  given  by 
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2  /  2 
x  X)  ■  ^  |110ux  -  [I20“xx  *  ’llV1  *  7  (^)  '>30\xx  ♦  2I21u, 
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6.2.3  Far  Field 


As  usual  Jn  denotes  the  unit  normal  vector  to  a  projected  flat  panel 
and  ig  and  j^  are,  respectively,  unit  vectors  along  the  x-  and  y-axis 
of  the  panel  coordinate  system  (equations  (7.2.10)  of  reference  1).  Define 
rQ  as  the  vector  from  the  origin  of  panel  coordinates  to  the  point  where 
velocity  is  being  evaluated  and  rQ  as  its  magnitude.  Certain  auxiliary 
vectors  are  needed: 
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The  far-field  expressions  for  the  source  velocities  are 

T/ (0 )  _  t2_I00  rQ 
2  r 
ro  0 


(77) 
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v  ^  20  U 
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w(ly)  _  tf  T  !o_t^  rT  t  ,  t  . 
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(79) 


The  fact  that  three  auxiliary  vectors  (76)  are  required  means  that  in 
effect  a  transformation  into  panel  coordinates  has  been  performed.  It  appears 
to  be  somewhat  faster  computationally  to  use  the  vector  far-field  forms  above 
rather  than  simply  truncate  the  intermediate  field  formulas  in  panel  coordi¬ 
nates.  However,  the  advantage  of  the  far-field  formulas  is  reduced  relative 

■f 

to  that  of  reference  1,  where  only  the  vector  D  is  required. 


6.3  Vorticity  Velocity  Formulas 

As  defined  in  equation  (59)  the  velocity  due  to  a  vorticity  distribution 
on  a  panel  is  expressed  as  a  sum  of  three  terms  which  are  denoted  by  super¬ 
scripts  0,  1,  and  2.  The  term  with  superscript  0  is  the  flat-panel 
constant-vorticity  contribution.  The  other  two  terms  account  for  surface 
curvature  and  vorticity  derivative  effects.  Specific  formulas  are  given 
below  for  these  velocities  in  terms  of  surface  curvatures  P,  Q  and  R  and 
derivatives  of  the  equivalent  dipole  distribution  u  (section  5.3).  Cal¬ 
culation  of  these  quantities  are  described  in  subsequent  sections.  It  is 
more  convenient  to  group  the  formulas  for  each  term  of  the  velocity  together 


than  to  organize  the  formulas  by  Near-,  Intermediate-,  and  Far-Fields,  as  is 
done  for  the  source  effects.  As  mentioned  above,  it  is  assumed  that  the 
formulas  below  are  being  evaluated  in  the  panel  coordinate  system. 

6.3.1  The  Superscript  tero  Velocity  Components 

It  turns  out  that  the  super  0  components  can  be  expressed  in  terms  of 
the  flat-panel  constant-source  velocities  (equation  (61))  as  follows 

=  -uxVz  (source) 

VSJ  =  -wyVz  (source)  (80) 

=  uxVx  (source)  -  uyVy  (source) 

These  are  evaluated  in  the  Near-  and  Intermediate-Fields  by  the  formulas  of 
reference  1  with  the  addition  of  first  area  moment  terms  in  the  Intermediate- 
Field  as  given  in  equation  (70).  In  the  Far-Field  the  first  terms  only  of 
the  Intermediate-Field  formulas  are  retained,  and  this  calculation  is  performed 
in  panel  coordinates. 

6.3.2  The  Superscript  One  Velocity  Components 

These  components  are  as  follows  in  the  Near-Field: 

Vix]  =  '^xx'ho  +  Voi}  +  MxH  +  2[KQ  ~  wyP)(‘Jll  +  xVy  (source)) 

+  (wxR  -  vyQ) (~ J02  +  ^vy  (source))] 

Viy^  =  '2^xy^l0  +  “yyW  +  vyH  +  2(("^XQ  +  yyP^‘J20  +  xVx  (source^ 

+  (-uxR  +  vyQ)(-J11  +  yVx  (source))]  (81) 

=  2^xx("J20  +  xVx  (source))  +  ^yy("J02  +  yVy  (source)) 

+  ^(-2^  +  yVx  (source)  +  xVy(source))] 

where 
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(82) 


H  =  P[J2Q  -  2 xVx  (source)]  +  2Q[J11  -  xVy  (source)  -  yVx  (source)] 

+  R[Jq2  -  2 yVy  (source)]  +  h 
and 

h  =  j  [Px2  +  2Qxy  +  Ry2] Vz  (source)  (83) 

For  small  z,  i.e.,  for  points  approaching  the  plane  of  the  projected 
flat  panel,  there  is  a  possible  singularity  in  equation  (83)  due  to  the  factor 
1/Z.  To  avoid  numerical  difficulty  a  "small -z"  test  is  employed  and  special 
procedures  employed,  if  z  is  small.  This  is  also  required  for  the  source 
formulas  used  in  both  the  base  method  and  the  present  method.  The  test  and 
the  procedures  are  described  in  section  7.8  of  reference  1.  In  particular 
Vz  (source)  is  set  equal  to  zero  if  the  projection  (x,y)  of  the  point  into 
the  plane  of  the  panel  lies  outside  the  panel. 


The  procedure  used  for  equation  (83)  uses  the  same  small -z  test  as 
reference  1,  and  if  it  is  satisfied  sets 


h  =  [Px‘  +  2Qxy  +  Ry^]-^  ,  (x,y)  outside  (84a) 

=  0  ,  (x,y)  inside  (84b) 


where  the  derivative  of  4>qq  is  given  in  equation  (7.7.9)  of  reference  1. 

The  reasoning  that  led  to  this  choice  is  as  follows.  Since  Vz  -*■  0  if  (x,y) 

is  outside,  it  is  legitimate  to  use  1 'Hospital's  rule  on  (83),  and  since 

Vz  =  ^OO’  e<luat'ion  (84a)  follows.  For  (x,y)  inside  the  panel  the  limit 

does  not  exist  for  arbitrary  values  of  x  and  y.  However,  the  assumption 

is  made  that  if  the  point  (x,y,z)  approaches  the  panel  it  approaches  along 

a  line  through  the  origin  that  does  not  lie  in  the  xy-plane.  Then  the  square 

o 

bracket  of  (83)  is  0 (z  )  and  (84b)  follows.  This  assumption  was  made  in  the 
derivations  of  section  5.0.  Furthermore,  it  is  physically  reasonable  that  if 
a  point  approaches  a  panel  it  does  so  at  the  panel  control  point  (i.e.,  the 
origin),  which  is  the  only  point  where  the  normal -velocity  boundary  condition 
has  been  applied. 


The  formulas  for  the  velocity  components  in  the  Intermediate  Field  are 
expressed  in  terms  of  the  direction  cosines  a,B,y  which  are  defined  in 
(7.6.2)  of  reference  1.  In  terms  of  the  auxiliary  potentials. 
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the  components  in  the  Intermediate  Field  are 
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V0)Z  =  2^xx^x<^10  ^20  ^  +  pxy{(x4>01  *11*  +  {y*10  *11 )}  +  pyy(y*01  *02)] 

where 


H  =  P*$Q  +  2Q*ft  +  R**2  (87) 

The  Far-Field  calculation  is  obtained  from  (85)  and  (86)  by  truncating 
(85)  with  the  first  term  containing  even-ordered  moments,  i.e.,  the  first  term 
if  p+q  is  even,  and  the  second  term  if  p+q  is  odd. 

6.3.3  The  Superscript  Two  Velocity  Components 

The  superscript  2  velocity  components  are  by  far  the  most  complicated 
analytically  in  the  near  field  and  a  considerable  number  of  secondary  quanti¬ 
ties  must  be  defined.  The  auxiliary  potentials  *mn  are  defined  by 

A  rf 

In  terms  of  these  potentials  the  near-field  formulas  are 
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(88) 


(89) 


38 


and 


V<2)  =  *  |P,JX^30  +  (2(\  +  yyP^21  +  (Ryx  +  2^uy^l2  +  R|Jy*03  +  2^xP  +  y<^px*20 
+  2[(xQ  +  yR)vx  +  (xP  +  yQ)wy]\|>11  +  2(xQ  +  yR)vy'l'02 
+  (Px2  +  2Qxy  +  Ry2)  (px^10  +  ^Ol^j 
To  evaluate  the  ^  the  following  quantities  are  defined 
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In  terms  of  these,  the  <|>  are 
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*21  =  Z  {  f32  lm32(q32  +  E32)V132)  +  q32932Vo32)]  +  f41  lm41(q41  +  E41)V141) 

*  W$411>  -77% L<32>  *  77% L<41’  *  L°2’  -  l(M)!  (104) 


p30 


=  z  - 


J_  2 

f32  ^m32E32 


—  a2  )V^32^  +  e  (E 
q32 '  V1  e32't32 


+  z2)Vp2)] 


+  f4l  C(m41E41  q41)Vl41)  +  e41(E41  +  z2)Vq41)] 

(41)) 


2  +  3m! 


'32  ,  (32) 


2  +  31^ 


(f32) 


m 


(f4l) 


372 


(105) 


A  "small-z"  test  is  required  for  ^00  from  (96).  From  (88)  and  (89)  it 


can  be  seen  that  t 


occurs  only  when  multiplied  by  a  factor  of 

iIiqq  is  not  really  singular,  and 


2  2  00 

z(Px  +  2Qxy  +  Ry  ).  Thus  the  1/z  term  of 
the  last  term  of  <|,qq  enters  only  in  the  combination  (83)  that  gives  h. 
Accordingly,  the  small-z  procedure  described  in  section  6.3.2  is  used  here 
as  well. 


Fortunately,  the  Intermediate- Field  formulas  for  these  components  are 
quite  simple.  The  auxiliary  potentials  emn  are  defined  as  suitably  trun¬ 
cated  expansions  of  the  integrals 


d£do 


Specifically, 


+  14aeIm+l,n+l 


+  B*m,n+l| 

+  (7e2-')Im,n«lJ  for  m  +  n  *  2 


(106a) 
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for  m  +  n  =  3 


(106b) 


V(dR  =  3z(Pe20  +  2^611  +  Re02^  (108) 

and 

v!z2)  =  (xux  +  yuy)V„R  -  32[pNe30  +  <2«\  +  Pl,y)e21  *  <R|1x  +  2Quy)612 

+  %«03'  (109) 

The  Far-Field  expressions  are  from  equations  (106)  through  (109)  by 
truncating  (106a)  with  the  first  term. 
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7.0  CALCULATION  OF  GEOMETRIC  QUANTITIES  FOR  A  PANEL 


7.1  General  Description  of  the  Surface- Curvature  Calculation 

As  illustrated  in  figure  1,  the  input  to  the  program  consists  of 
coordinates  of  a  set  of  points  that  define  the  body  surface.  The  order  of 
input  is  such  that  successive  input  points  define  a  curve  in  the  body  surface 
called  an  N-line.  When  such  a  curve  has  been  completely  defined,  the  input 
logic  notes  the  fact  and  subsequent  points  define  an  adjacent  N-line.  Four- 
cornered  panels  are  formed  from  points  on  two  adjacent  N- lines  and  thus  there 
is  a  "strip"  of  panels  between  adjacent  N-lines  (figure  1).  Thus  the  ordering 
of  surface  panels  is:  all  panels  of  one  strip  from  first  to  last  followed  by 
all  panels  of  the  next  strip,  etc.  This  ordering  of  the  panels  is  important 
both  for  the  curvature  calculation  of  this  section  and  the  source  derivative 
calculation  of  section  8.1. 

Suppose  the  geometric  quantities  for  a  typical  panel  are  being  calculated 
(the  "panel  in  question"  in  figure  4).  A  least-square  plane  is  computed  for 
the  four  input  points  that  form  the  corners  of  this  panel  in  the  manner  des¬ 
cribed  in  references  1  and  2.  A  coordinate  system,  £,n,t»  is  constructed 
in  this  plane  having  the  c-axis  normal  to  the  plane  and  having  as  origin  the 
point  whose  coordinates  are  the  averages  of  those  of  the  four  input  points. 

The  four  input  points  are  transformed  into  the  S,o,C-coordinate  system.  Also 
transformed  into  this  system  are  the  additional  input  points  that  form  the 
corners  of  the  four  panels  adjacent  to  the  one  in  question.  As  shown  in 
figure  4,  each  adjacent  panel  provides  two  additional  input  points,  because 
the  other  two  are  identical  with  two  of  those  forming  the  corners  of  the  panel 
in  question.  The  expression  defining  the  paraboloidal  panel  is  assumed  in 
the  form 

C  =  ZQ  +  Ac  +  Bn  +  R<i2  +  2Qcn  +  Rn2  (110) 

The  six  parameters  zQ,  A,  B,  P,  Q,  and  R  are  determined  from:  (a)  requir¬ 
ing  expression  (110)  to  pass  through  the  transformed  corner  points  of  the 
panel  (four  conditions),  and  (b)  requiring  expression  (110)  to  agree  as  well 
as  possible  in  a  least  squares  sense  with  the  eight  additional  input  points 
for  the  adjacent  panels.  This  last  is  performed  with  two  independent  param¬ 
eters  and  thus  represents  two  additional  conditions  -  a  total  of  six  as 
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required.  If  the  panel  in  question  is  the  first  or  last  of  a  strip,  there 
are  only  three  adjacent  panels.  This  is  also  true  for  most  panels  on  the  first 
and  last  strips  of  a  body.  The  first  and  last  panels  of  the  first  and  last 
strips  have  only  two  adjacent  panels.  Thus,  there  may  be  only  six  or  four 
additional  input  points  and  the  least-squares  procedures  are  used  accordingly. 
There  is  no  difficulty  since  only  two  parameters  are  being  adjusted.  As 
described  in  the  next  section,  it  is  necessary  to  perform  this  calculation  in 
double-precision  arithmetic. 

It  should  be  noted  that  the  form  (110)  is  not  equivalent  to  the  form  (6) 
required  by  the  derivations  of  section  5.0.  This  is  because  the  least-square 
plane  is  not  the  plane  of  the  flat  projected  panel.  Rather  the  least-square 
plane  is  a  temporary  device  needed  to  perform  calculations.  Once  the  form 
(110)  has  been  obtained,  the  constant  and  linear  terms  are  eliminated  by 
translation  and  rotation  of  coordinates.  The  new  plane  obtained  by  these 
translations  and  rotations  is  that  of  the  flat  projected  panel. 


7.2  Specific  Formulas  for  the  Surface-Curvature  Calculation 


The  initial  portion  of  this  calculation  is  identical  with  that  of  the 
flat-panel  procedure  of  reference  1.  Specifically,  a  "chord-plane"  flat  panel 
is  constructed  having  the  average  point  as  origin  and  having  the  c  axis  in 
the  normal  direction  as  defined  in  reference  1.  Equations  (7.2.1)  through 
(7.2.15)  of  reference  1  are  carried  out.  Then  the  four  input  points  used  to 
form  the  panel  are  transformed  into  panel  coordinates  by  the  usual  equations. 
Let  nK>  K  =  1,  2,  3,  4  denote  the  panel  coordinates  of  the  input 
points.  The  numbering  is  counterclockwise  around  the  panel  as  viewed  from  the 
negative  ;-axis  (figure  20  of  reference  1).  The  following  quantities  are 
now  defined. 


d  =  2  ^3  ^4^ 


e  2  (?2  “  *1 ) 


(e,  ♦  -  e3  -  s4) 


h 


(111) 


e  “  T  (”1  “  *2  +  n3  ”  n4^ 
f  ■  }  t'l  -  {2  *  '3  "  '«> 
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Then  necessary  derived  quantities  are 


ZP  "  "7  (d  +  e  +  2q ' 

Zg  =  -e(d  -  e)  -  2qh 
zR  =  -(h2  +  e2) 


h(d  +  e) 


n  -  d  -  e 
yA  '  2h(d  +  e) 

%  =  h(d  +  e) 


«p-I 
,  _  1 


AA  d  +  e 

r  -  e  ~  d 
CAB  "  r+T e 

n  =  f  e  ~  d 

uo  T  rnr 


2de  -  *»ctdh-  e> 


(112) 


DR  =  2he 
CBA  ■  2o + 

CBB  *  I  <"2  -  '2! 

2o--rf 

Ep  =  d2  —  e2 


A  2x2  set  of  linear  equations  is  solved  for  three  different  right 
sides  to  obtain  the  quantities  AQ>  BQ,  AR,  BR,  Ap,  Bp,  namely 


CAAAo  +  CABBo  =  Do 
CBAAo  +  CBBBo  =  Eo 


(113a) 


caaar  +  cabbr  =  °R 
cbaar  +  cbbbr  =  0 


(113b) 
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(113c) 


caaap  +  cabbp  =  0 
cbaap  +  cbbbp  =  EP 


Now  if  the  quantities  A,  B,  Q  and  zQ  in  (110)  are  taken  in  the  form 


A  =  A0  +  ApP  +  ArR 
B  =  Bq  +  BpP  +  BrR 
Q  =  Q0  +  Qaa  +  QbB  +  QpP 


(114) 


ZPP 


+  zqQ  +  zrr 


then  expression  (110)  passes  through  the  input  points  f;^,  nK»cK,  K  =  1,2, 3,4 
for  arbitrary  values  of  P  and  R.  These  two  quantities  are  adjusted  to  give 
a  least-square  fit  to  the  surrounding  panels. 


If  the  quantities  A,  B,  Q  and  zQ  are  replaced  in  (110)  by  their 
representations  in  (114),  the  result  is 

C  =  D(c,n)  +  EU,n)P  +  F(e,n)R  (115) 

where 

D(C,n)  =  zq(Q0  +  QaAo  +  QbBq)  +  AqC  +  B0n  +  (2Qq  ♦  2Q^0  +  2QBBQ)5n 

E(5,n)  =  zp  +  zQ(QAAp  +  QgBp  +  Qp)  +  Ape  +  Bp^  +  52  +  (2QAAp  +  2QgBp  +  2Qpkn 

2  (116) 

F(c,n)  =  zq(Qaar  +  QgBR)  +  zR  +  Ar?  +  BRn  +  n  +  (2QAAR  +  2QgBR)Cn 

Figure  4  shows  the  four  panels  surrounding  the  "panel  in  question"  whose 
geometry  is  being  calculated.  For  the  interior  panel  shown,  the  four  surrounding 
panels  are  defined  by  a  total  of  eight  input  points  that  are  not  also  input 
points  defining  the  panel  in  question.  As  detailed  in  section  7.1,  panels  at  the 
edge  of  a  section  may  have  2  or  3  adjacent  panels,  which  correspondingly  may 
have  4  or  6  input  points.  The  input  points  of  the  surrounding  panels  are  trans¬ 
formed  into  the  coordinate  system  of  the  panel  in  question.  The  resulting  coordi¬ 
nates  are  £m»  nm»  cm.  m  =  1 ,  ....  8  (or  6  or  4).  The  least-square  procedure  is 
accomplished  by  defining  the  quantities 
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I 


( 

1  i 


8 

CRR 

-E 

IF^m 

•\»2 

m=l 

8 

- 

CPP 

■E 

‘E<V 

’nm)]2 

m=l 

8 

CRP 

=  CPR 

II 

M 

^m’V^m’V 

m=l 

8 


L1  D^m*nm^F^m*nm^ 

m=l 

8 

L2  “  D^m*nm^E^m’T1m^ 

m=l 

after  which  P  and  R  are  solved  for  from  the  2x2  system  of  equations 


CRRR  +  CRPP  =  L1 
cprR  +  CppP  =  l2 


018) 


If  the  panel  in  question  is  on  the  edge  of  a  section,  the  summations  in 
(117)  extend  to  4  or  6  as  appropriate  rather  than  to  8.  It  was  found  that 
on  wings  where  panel  aspect  ratios  of  50  are  not  uncommon,  the  fact  that 
four  of  the  points  Um,nm,tm)  have  much  larger  £n  coordinates  than  the 
other  four  can  lead  to  ill-conditioning  of  equations  (118)  due  to  round-off 
error.  Accordingly,  equations  (117)  and  (118)  are  calculated  using  double¬ 
precision  arithmetic. 


Once  P  and  R  are  known,  all  other  parameters  of  (110)  can  be 
obtained  from  (114). 


If  one  side  of  the  panel  is  much  smaller  than  the  other,  some  instability 
can  result  in  equations  (111),  (112),  and  (113).  Accordingly,  in  this 
case  special  "zero  length  side"  formulas  are  used: 

f  =  e  =  0 

A0  =  Ar  =  Ap  =  0  (119) 

Bo  =  BR  =  0 


L 
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7.3  General  Description  of  the  Calculation  of  the  Control  Point  and  the 

Projected  Panel 

All  of  the  calculations  of  sections  7.1  and  7.2  are  carried  out  in  a 
coordinate  system  whose  Sn-plane  is  the  "chord  plane",  which  is  also  the 
plane  of  the  flat  panel  of  reference  1.  However,  this  is  not  the  plane  of 
the  projected  flat  panel  of  the  present  method  as  described  in  section  5.0. 
Among  other  things  the  projected  flat  panel  is  tangent  to  the  true  curved 
surface.  The  tangency  point  is  both  the  control  point  of  the  panel  and  the 
origin  of  the  panel  coordinate  system.  It  remains  to  compute  the  tangency 
point,  the  normal  direction  at  the  tangency  point,  and  the  geometric  quanti¬ 
ties  associated  with  the  flat  projected  panel.  These  quantities  may  be 
defined  in  several  ways,  all  of  which  are  equivalent  to  the  order  of  accuracy 
of  the  overall  calculation.  The  selections  made  for  the  present  method  are 
outlined  below. 

The  control  point  is  taken  as  that  point  of  the  surface  (110)  that  is 
nearest  the  origin,  and  the  normal  direction  is  that  from  the  origin  to  this 
nearest  point.  Thus,  to  the  order  of  accuracy  being  considered,  the  control 
point  is  the  point  on  the  true  curved  surface  nearest  the  average  point  of  the 
four  input  corner  points.  The  size  and  shape  of  the  projected  flat  panel  are 
assumed  identical  to  those  of  the  "chord  plane"  flat  panel  calculated  by  the 
procedure  of  reference  1.  Thus  the  flat  projected  panel  is  determined  solely 
by  the  coordinates  of  the  control  point  and  the  unit  vectors  defining  the  axes 
of  the  panel  coordinate  system,  whose  c-axis  defines  the  normal  direction. 


7.4  Specific  Formulas  for  the  Calculation  of  the  Control  Point  and  the 
Projected  Panel 

The  derivatives  of  expression  (110)  are 


?£  =  A  +  2(PC  +  Qn) 
Cn  =  B  +  2(Q5  +  Rn) 


(120) 


The  point  on  the  surface  nearest  the  origin  is  characterized  by  the  condition 
that  the  vector  from  the  origin  is  parallel  to  the  local  normal  vector,  or 

Ul  +  n j  +  cU»n)£]  x  (-c^i  -  c^j  +  k]  =  0  (121) 
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This  last  is  equivalent  to  the  two  scalar  equations 


G(?,n)  =  £  +  ?U»n)CgU,n)  55  0 
H(c,n)  =  n  +  ?(£»n)cr)(C»n)  =  0 


(122) 


The  nonlinear  equations  (122)  are  solved  by  Newton-Raphson  iteration.  The 
derivatives  are 

G?  =  1  +  +  2P c 

Gn  =  H?  =  +  2Q?  0  23) 

H  =  1  +  £  +  2RC 
n  n 

Let  5k»nk  denote  the  approximations  to  the  solutions  of  (122)  obtained 
in  the  k-th  iteration.  The  next  approximation  is  obtained  by  solving  the  2x2 
set  of  linear  equations 


G5(ek,\)(5k+1  +  Gn(5k’nk)(nk+1  =  '^k.V 

Vek,nk^5k+1  ~  Sk)  +  Hn(^|c»rik)(nk+1  -  nk)  =  ’HUk>ok) 


(124) 


Initial  values  are  taken  as  =  n0  =  0.  The  iteration  converges  very 
rapidly  to  the  nearest  point  with  coordinates  The  local  normal 

vector  at  this  point  (which,  of  course,  is  parallel  to  the  vector  from  the 
origin)  has  components  y-|.y2*Y3  in  this  coordinate  system,  where 


*1  =  '^(S^nJ/Y 

y2  =  'cn(c»  ’0/y  025) 

_  S  2  2 

Y3  ~  -  Yt  -  Y2 

Y  =Vl "  +  \\  +  ^  (126) 


The  above  vector  is  the  unit  normal  vector  to  the  projected  flat  panel  and 
defines  the  third  axis  of  its  coordinate  system.  It  is  now  necessary  to 
construct  unit  vectors  parallel  to  the  other  two  axes  of  this  system.  Pre¬ 
liminary  quantities  are 


Y1 


(127) 
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If  both  y-|  and  y2  are  essentially  zero,  set  y^  =  y^  =  f7/2.  The  compon¬ 
ents  of  the  vector  parallel  to  the  first  axis  of  the  projected  flat  panel's 
coordinate  system  are  a-1.a2.a3  where 

al  =  ^Y2^2  +  y3^1^ 

a2  =  ^y3  ~  1 )yiy2  028) 

a3  -  -Y]  VI  -  y3 


Finally,  the  components  e-|  6^3  of  the  unit  vector  parallel  to  the  second 
coordinate  axis  are  obtained  by  taking  the  cross  product  of  the  vector  (125) 
and  the  vector  (128).  The  matrix 


would  thus  be  suitable  for  transforming  vectors  and  points  between  the 
coordinate  systems  of  the  projected  flat  panel  and  the  "chord-plane"  panel. 
As  defined  in  equation  (7.2.10)  of  reference  1,  the  transformation  matrix  of 
the  "chord-plane"  panel  is. 


~all 

a12 

a13~ 

la(old)]  = 

a21 

a22 

a23 

-a31 

a32 

a33  - 

O30) 


Thus  the  transformation  matrix  for  the  projected  flat  panel  is  obtained  as 
the  matrix  product 


[a(new)]  =  [A]  [a(old))  (131) 

The  rows  of  matrix  (131)  are  unit  vectors  parallel  to  the  axes  of  the 
coordinate  system  of  the  projected  flat  panel  expressed  in  the  reference 
coordinate  system.  In  particular,  the  third  row  is  the  normal  vector.  This 
matrix  is  used  in  the  manner  of  reference  1  to  transform  points  and  vectors. 
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The  reference  coordinates  of  the  control  point  are  obtained  from 


xo  =  xav  +  allc»  +  a21n»  +  a31c® 

yo  =  yav  +  a12^»  +  a22n®  +  a32?»  (132) 

2o  =  zav  +  a13^»  +  a23n®  +  a33<® 

The  quantities  on  the  right  side  of  (132)  are  those  of  (130)  not  (131). 

However,  this  is  the  last  use  of  matrix  (130)  which  is  now  discarded  in  favor 
of  (131)  just  as  the  average  point  (*av»  yav»  *av)  is  discarded  in  favor  of 
(xQ,  yQ,  zQ).  Because  of  the  way  the  projected  flat  panel  has  been  defined, 
these  are  the  only  required  replacements.  All  other  quantities  that  have  been 
computed  for  the  "chord-plane"  panel  are  valid  also  for  the  projected  flat 
panel,  specifically  the  slopes  and  corner-point  coordinates,  equations  (7.2.12) 
and  (7.2.15)  of  reference  1.  The  remaining  geometric  quantities  of  the  pro¬ 
jected  flat  panel  are  calculated  in  the  manner  of  reference  1  beginning  with 
equation  (7.2.20). 

It  could  be  argued  that  the  entire  process  of  section  7.0  could  be  iter¬ 
ated.  Using  the  new  panel  coordinate  system,  the  procedure  of  sections  7.1 
and  7.2  would  give  somewhat  different  values  of  the  coefficients  in  (110) 
which  could  in  turn  be  used  in  the  calculation  of  the  present  section.  How¬ 
ever,  if  panel  dimensions  are  small  compared  to  physical  dimensions  of  the 
body,  as  is  assumed  throughout  the  present  method,  s  is  small  everywhere 
on  the  panel  and  the  rotations  implied  by  matrix  (129)  are  also  small.  It 
can  be  shown  that  the  changes  to  be  realized  by  iteration  are  of  higher  (third) 
order  in  panel  dimensions,  and  it  is  consistent  with  the  other  approximations 
that  have  been  made  not  to  iterate.  This  course  has  been  adopted  for 
simplicity. 
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8.0  THE  SOURCE  DERIVATIVE  TERMS.  ASSEMBLY  OF  THE  MATRIX  OF 
INFLUENCE  COEFFICIENTS 


8.1  The  Least-Square  Procedure.  Geometric  Constants 

As  stated  in  section  6.2,  the  basic  source  velocity  formula  (60)  contains 

coefficients  a  and  a  ,  which  are  the  derivatives  of  the  source  density 
x  y 

with  respect  to  the  panel's  coordinate  directions.  It  is  not  intended  that 
these  be  additional  unknowns.  Instead,  they  are  expressed  in  terms  of  the 
unknown  values  of  source  density  at  the  control  points  of  the  surrounding 
panels.  Thus,  ultimately  the  values  of  source  density  at  the  control  points 
of  the  panels  are  the  only  unknowns.  Just  as  for  the  curvature  calculation 
of  section  7.1,  the  source-derivative  procedure  is  slightly  different  for  the 
first  and  last  panels  of  a  strip  and  for  panels  of  the  first  and  last  strips. 
However,  the  modifications  are  quite  straightforward.  In  the  initial  discus¬ 
sion  it  is  assumed  that  the  panel  on  which  source  derivatives  are  being 
evaluated  (the  panel  in  question)  has  adjacent  panels  on  all  four  sides  as 
shown  in  figure  5.  For  the  purposes  of  the  present  discussion  only,  the  con¬ 
trol  points  of  the  adjacent  panels  are  numbered  K  =  0,  1,  2,  3,  4  as  shown 
in  figure  5,  where  0  denotes  the  element  in  question.  These  control  points 
are  transformed  into  the  coordinate  system  of  the  panel  in  question.  Let  the 
£  and  n  coordinates  of  these  control  points  be  £qK,  r^,  K  =  0,  1,  ....  4 
and  the  values  of  source  density  at  the  control  points  be  a^,  K  =  0,  1 ,  . . . ,  4 
(evidently  £qq=  ngo=  0).  The  source  density  on  the  panel  in  question  is 
assumed  in  the  form 


a(£,n)  =  °0  +  ox£  +  ayn  (133) 

This  is  fitted  in  a  least-square  sense  to  the  values  of  source  density  at  the 
control  points  of  the  adjacent  panels.  That  is, the  source  density  variation 
on  the  panel  in  question  is  a  least-square  plane  through  neighboring  values. 
This  process  expresses  the  source  derivatives  on  the  panel  in  question  in  terms 
of  the  unknown  values  of  source  density  on  the  adjacent  panels  in  the  form 
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(134) 


£  4y,< 


K 


where  the  coefficients  cix^ 

K 


K=0 


and  are  purely  geometrical. 


Specifically,  let  M  be  the  number  of  control  points  involved  in  the 
least-square  fit*  Thus  M  =  5  for  an  interior  panel  as  shown  in  Figure  5, 
but  M  =  4  for  most  panels  of  the  first  and  last  strips  and  for  the  first 
and  last  panels  of  interior  strips.  For  the  first  and  last  panels  of  the  first 
and  last  strips  M  =  3.  The  following  quantities  are  defined 

M-l  M-l 


X  = 


’OK 


Y  = 


K=0 


'OK 


K=0 


M-l 


*•£ 
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K=0 
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-£ 
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•£ 

K=0 
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a  =  S 


r 

M 


h  -  t  XY 

b  ~  T~r 


c  =  U  - 


d  =  ac  —  b 


Then  the  desired  coefficients  are 


r(x)  _  c  _  _  b  ,  b  Y  cX 

K  d  ^OK  d  n0K  c  M  d  M 


(•(^  =  -  -  f  +  —  n 

Ct/  d  50K  a  n 


a 

d 


,  b  X  _  a  Y 
OK  d  M  dfi 


K  =  0,  1 . M 


(136) 


There  are  2M  "source  derivative  coefficients,"  a  total  of  10  for  an  interior 
panel.  These  are  calculated  once  and  for  all  and  stored  with  the  other 
geometric  quantities  defining  each  panel. 


8.2  Logic  of  the  Assembly  Procedure 

In  the  base  method  the  velocity  induced  by  a  panel  depends  only  on  the 
source  density  on  that  panel  and  thus  the  "influence  coefficients"  for  that 
panel  are  calculated  solely  from  that  panel's  geometry.  The  essentially  new 
feature  of  the  source  derivative  procedure  is  that  the  velocity  induced  by  a 
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panel  depends  on  the  value  of  source  density  at  the  control  point  of  that  panel 
and  also  on  the  values  of  source  density  at  the  control  points  of  adjacent 
elements.  Similarly,  the  velocities  induced  by  adjacent  elements  depend  on 
the  source  density  on  the  panel  in  question.  Thus  the  "influence  coefficients" 
for  a  panel  depend  not  only  on  the  geometry  of  that  panel  but  also  on  the  geom¬ 
etry  of  adjacent  panels  and  the  assembly  of  the  influence  coefficient  matrix  is 
more  complicated. 


Let  the  panels  be  numbered  consecutively  in  the  order  they  have  been 

formed.  Thus,  reference  is  made  to  the  i-th  panel  and  to  the  j-th  panel  where 

both  i  and  j  range  from  1  to  N.  Another  way  of  stating  the  essentially 

new  feature  above  is  that  a  distinction  must  be  made  between  the  effect  of  the 

j-th  panel  and  the  effect  of  the  j-th  value  of  source  density,  whereas  these 

two  effects  are  identical  in  the  base  method.  Let  V*.  be  the  velocity 

1 J  J 

induced  at  the  i-th  control  point  by  the  j-th  panel  and  V^.  be  the  velocity 
induced  at  that  point  by  the  j-th  value  of  source  density.  Then  in  the  nota¬ 
tion  of  section  6.2  and  the  present  section, 

)(0)  +  ^(P), 


7*  - 


1J 


=  [V' 


'P  +  2V(Q)Q  +  V(R)R  +  c(x)V(lx)  +  c(y)V(1y)]a 

O  0  0 


M-l 


T,  t4x)v(lx)  +  4y^ly)]oF 


(137) 


Notice  that  subscripts 
(137)  for  simplicity, 
and  o-i,  o 2»  Og,  and 


K=1 


i  and  j  are  omitted  on  the  right  side  of  equation 

)  is  o. 
in  (137) 


In  the  overall  numbering  scheme,  aQ  in  (137)  is  a 


have  subscripts  near  j.  All  velocities 

The  curvatures  P,  Q,  R  and 


1 »  “2*  u3*  u4 

depend  only  on  the  geometry  of  the  j-th  panel 
^  and  cj^ 

calculated  they  can  be  associated  with  the  j-th  panel  only. 


the  coefficients  ciA/  and  c„J 1  depend  on  the  surrounding  panels,  but  ^nc? 


Consider  now  the  i-th  row  of  the  matrix  V. .,  which  expresses  the 

•  J 

effects  of  the  various  values  of  source  density  at  the  i-th  control  point. 

The  first  bracketed  term  in  (137)  is  an  effect  of  ai  and  is  added  to  the 
j-th  location  of  the  row.  The  four  terms  in  the  sunmation  of  (137)  represent 
effects  of  other  values  of  o  and  must  be  added  to  other  locations.  Refer¬ 
ring  to  figure  5,  it  can  be  seen  that  the  panels  numbered  1  and  2  are  on  the 
same  strip  as  the  panel  in  question  and  thus  represent  effects  of  the  proceed¬ 
ing  and  succeeding  values  of  a.  In  particular,  value  1  is  associated  with 
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Oj _i  and  value  2  with  and  the  relevant  terms  of  (137)  are  added  to 

those  locations.  Panels  3  and  4,  however,  are  on  adjacent  strips.  Suppose 
there  are  E  panels  on  each  strip.  Then  value  3  is  associated  with  o.  F 
and  value  4  with  Oj+£,  and  the  relevant  terms  of  (137)  are  added  to  these 
locations  of  the  row. 

If  the  matrix  of  influence  coefficients  was  formed  row  by  row,  or  if  it 
could  be  completely  contained  in  high-speed  storage,  the  above  procedure 
would  be  fairly  simple  (as  indeed  it  is  in  two  dimensions).  However,  the 
matrix  is  formed  column  by  column,  and  some  rather  complicated  logic  is 
required.  By  retaining  always  three  columns  of  the  matrix  in  high-speed  stor¬ 
age,  the  terms  in  (137)  that  belong  in  columns  j-1,  j,  and  j+1  can  be  put 
into  the  proper  locations  immediately  and  the  results  stored  on  the  main 
matrix  tape  (or  disk)  in  the  usual  way.  However,  the  j-E  effects  must  be 
stored  on  a  second  tape  (or  disk)  and  the  j+E  effects  on  a  third  tape. 

Thus,  when  the  usual  calculation  is  over  and  all  values  of  j  have  been  con¬ 
sidered,  the  result  is  three  tapes  each  containing  an  almost  full  matrix. 

(The  only  deviations  from  full  matrices  are  that  the  matrix  on  the  second 
tape  has  zeros  in  the  last  E  columns  and  that  on  the  third  tape  has  zeros 
in  the  first  E  columns.)  The  three  matrices  are  brought  into  high-speed 
storage  column  by  column,  and  corresponding  columns  are  added  location  by 
location  to  produce  the  final  V.,  matrix.  To  save  time  the  induced  normal 

•  *  J 

velocity  matrix  is  also  computed  at  this  time  by  taking  dot  products  with  the 
unit  normal  vectors  to  the  panels. 
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9.0  THE  GLOBAL  VORTICITY  VARIATION 


9.1  General  Description  of  the  Procedure 

The  formulas  of  section  6.3  give  velocity  components  at  a  point  due  to  a 
vorticity  distribution  on  a  curved  panel  in  terms  of  the  first  and  second 
derivatives  of  the  equivalent  dipole  distribution  u(c,n)  evaluated  at  the 
panel  "center,"  i.e.,  the  origin  of  panel  coordinates.  Thus  there  are  five 
parameters  per  panel,  ux,  uy,  uxx,  vyy,  and  these  permit  a  general 
quadratic  dipole  variation.  These  derivatives  are  determined  by  the  algo¬ 
rithm  adopted  for  global  vorticity  or  dipole  variation. 

As  mentioned  in  section  4.0,  vorticity  is  an  auxiliary  singularity  in 
the  present  method,  and  its  main  purpose  is  to  satisfy  the  Kutta  condition 
along  a  wing  trailing  edge.  Thus  a  global  variation  algorithm  is  chosen  that 
expresses  the  vorticity  distribution  over  the  entire  body  (or  the  equivalent 
dipole  distribution)  in  terms  of  a  number  of  parameters  equal  to  the  number 
of  trailing-edge  segments  at  which  the  Kutta  condition  is  applied  (figure  1). 

In  the  logical  scheme  of  the  present  method  there  are  three  aspects  to  the 
global  variation:  (1)  global  variation  along  an  N-line  (figure  1),  (2)  local 
variation  along  an  N-line  over  a  single  panel,  and  (3)  local  variation  between 
consecutive  N-lines  across  the  intervening  strip  of  panels.  The  global  vari¬ 
ation  across  the  N-lines  ("spanwise")  is  not  part  of  the  global  variation 
algorithm  but  is  determined  by  the  Kutta  condition. 

The  global  variation  along  an  N-line  is  the  most  basic  assumption  and 
the  one  that  permits  the  most  flexibility.  The  distribution  of  the  equivalent 
dipole  strength  is  taken  as  a  particular  monotonic  function  of  arc  length  along  the 
N-line,  beginning  at  the  trailing  edge,  coming  forward  to  the  leading  edge 
along  the  lower  surface,  proceeding  along  the  upper  surface  back  to  the 
trailing  edge,  and  continuing  into  the  wake.  Section  7.3  of  reference  1 
gives  a  very  complete  description,  which,  while  it  is  specialized  to  the  case 
of  a  particular  distribution  function,  illustrates  the  general  case  very  well. 

(As  discussed  there  the  N-line  may  be  traversed  in  the  opposite  sense  also, 
i.e.,  upper  surface  first.)  It  is  understood  that  the  actual  variation  of 
equivalent  dipole  strength  along  the  N-line  is  the  product  of  the  distribution 
function  and  an  initially  unknown  multiplicative  constant  whose  value  is  a 
measure  of  the  local  "bound"  vorticity  strength.  Values  of  the  multiplicative 


constants  on  all  N-lines  are  obtained  as  solutions  of  the  simultaneous  equa¬ 
tions  that  express  the  Kutta  condition  at  the  trailing-edge  segments  of  the 
strips  (reference  1  and  section  10D).  In  the  usual  case  the  equivalent  dipole 
strength  on  the  wake  portion  of  the  N-line  is  taken  as  constant  and  equal  to 
its  value  at  the  upper  surface  trailing  edge.  As  described  in  the  paper  form 
of  reference  1,  the  base  method  has  two  options  for  the  global  dipole  varia¬ 
tion  along  an  N-line:  linear  and  cubic,  which  correspond,  respectively,  to 
a  constant  and  a  parabolic  "chordwise"  variation  of  bound  vorticity. 

Contradictory  as  it  may  first  appear,  the  assumptions  of  the  global 
variation  algorithm  also  determine  the  local  variation  over  a  panel  of  the 
vorticity  or  dipole  strength.  Since  the  formulas  of  section  6.3  permit  a 
linear  vorticity  variation  (quadratic  dipole  variation),  no  higher-order 
variation  can  be  considered.  On  the  other  hand,  that  order  of  variation 
should  be  used  to  obtain  the  full  power  of  the  present  method.  In  the  base 
method  of  reference  1,  the  dipole  variation  over  the  "chord"  of  the  panel  is 
linear,  while  its  "spanwise" variation  may  be  either  constant, linear, or 
quadratic. 

Because  only  a  limited  amount  of  effort  could  be  devoted  to  a  preliminary 
version  of  the  present  method,  it  was  decided  to  employ  the  global  variation 
algorithm  of  the  base  method  of  reference  1  (including  the  refinements  of  the 
paper  version),  even  though  some  of  the  capability  would  not  be  utilized. 

It  Is  intended  to  develop  a  more  sophisticated  algorithm  at  some  future  time. 

9.2  Formulas  for  Constant  Chordwise  Vorticity 

In  this  option,  the  equivalent  dipole  strength  p  is  assumed  to  vary 
linearly  along  an  N-line  in  the  form 

p  =  Bs  (138) 

where  B  is  the  unknown  multiplicative  parameter  and  s  is  arc  length 
measured  from  the  trailing  edge.  Over  an  individual  panel  the  variation  is 

P  =  A  +  B5  -  B(h  +  £)  (139) 

where  £  has  its  usual  meaning  as  a  panel  coordinate  (figure  20  of  reference  1) 
and  h  is  the  total  arc  length  along  the  N-line  up  to  the  n-axis  (c  =  0) 
of  panel  coordinates.  The  two  N-lines  bounding  the  panel  are  denoted  the 
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"first"  and  the  "second"  and  quantities  associated  with  them  are  identified 
by  subscripts  F  and  S,  respectively.  Then  the  dipole  strength  on  the  panel 
is  (compare  equation  (7.3.7)  of  reference  1) 


_  BF  BSr  .  BFftF  BShS  .  BSnl  BFn3  _  .  BShSnl  BFhFn3 
- —  ^  + - w - n  + - w - ?  + - w - 


+  c(Bp  -  Bs)(n  -  n3)(n  -  n-j )  +  d(Bp  -  Bs)f;< 


(140) 

where  w  is  the  panel  width  and  n>|  and  n3  are,  respectively,  the 

n-coordi nates  of  the  first  and  second  N-lines.  By  comparison  with  equation  (50) 

of  section  5.3  the  above  gives  the  required  five  parameters  as 

BSnl  ~  BFn3 
px  = - w - 


_ 


Bphp  —  Bjh^ 


w 


c(Bp  B^)(ti^  +  n3 ) 


(141) 


^xx  =  d(BF  "  BS) 

.  1  BF  ~  BS 
pxy  ”  2  w 

Pyy  =  C(Bp  -  BS) 


For  future  applications  provision  has  been  made  for  the  constants  c  and  d 
for  each  panel,  but  in  .the  current  form  of  the  present  method  d  is  always 
zero  and  c  is  nonzero  only  on  wake  panels  as  explained  in  reference  1. 


From  (140)  it  is  clear  that  the  equivalent  dipole  distribution  has  the 

form 


y  =  Bppp  +  B$ps 


(142) 


where  Up  and  us  are  independent  of  the  values  of  B  and  can  be  obtained 
directly  from  (140).  In  particular  up  is  a  dipole  distribution  that  is 
zero  on  the  second  N-line  and  has  a  unit  derivative  on  the  first.  The 
velocity  due  to  the  vorticity  distribution  on  the  panel  also  can  be  written 

V„  -  ♦  Vu(S)Us  <’«> 

where  the  velocities  Vu(F)  and  Vu(S)  correspond,  respectively,  to  equivalent 
dipole  distributions  pp  and  p^.  As  in  equation  (59)  of  section  5.3,  each 
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of  these  two  velocities  is  the  sum  of  three  terms 


V  (F)  =  V^(F)  +  V^V)  +  V*2V) 

(i)  tol  W 

V  (S)  ■  ?°)(S)  +  V(1)(S)  +  V'2)(S) 

U)  (*)'  <jJ  '  '  fa)'/ 


(144) 


Specific  assembly  formulas  for  the  components  of  the  velocities  in  (144) 
are  given  below. 


9.2.1  Near-Field  Formulas  for  Panels  on  the  Body 
The  superscript  zero  components  are 


Vix}(S)  “  -_w  Vz  (source) 


"1 


v2)(p)  *  1?  Vz  (source) 

hs 

w  c(nl  +  T13) 
hF  I 

^-~C(T11  +  'A 


uX 
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Vz  (source) 


V2  (source) 


(145) 


vi2>tS)  =  5T  Vx  <s°urce)  -  IS1-  «("!  +  "3> 

vlz1(F)  ■  -?vx  (souto)  *  1 +  y 


(source) 


Vy  (source) 


where  all  quantities  on  the  right  have  their  usual  meanings.  Near-field 
formulas  are  (7.7.7)  of  reference  1. 

Define  auxiliary  quantities  by 

J1  =  "J20  +  xVx  (source) 

JS  =  -J,,  +  xVw  (source) 

d  11  y  (*<46a ) 

J3  =  'Jll  +  yVx  (source) 

J4  =  "J02  +  yVy  (source) 

The  Jmn  have  the  meanings  appropriate  to  the  source  formulas.  Specifically, 
J2Q  is  given  by  equation  (64),  by  (7.7.12)  of  reference  1,  and  JQ2  by 


J02  =  H02  “  zVz  (source) 


(146b) 


where  Hgg  is  given  by  (7.7.15)  of  reference  1. 

In  terms  of  these,  the  superscript  one  components  are 

Hi  1  Oi  rii 

vix  (s)  =  2<d*in  +  SI  +  4  H  +  2  4  «WS  ♦  RJJ) 
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fhF  1 

+  [jr-c^l  +  n3)J  (H  +  2PJf  +  2QJ*) 


V!z)<S>  *  -2[dJi  +  CJ4  ♦  b  (J2  *  J3» 
V«<F)  -  2tdJf  +  cj*  ♦  ij  (0|  *  J*)l 

where  H  is  defined  by  (82). 


In  the  evaluation  of  the  superscript,  two  components,  V^R  i: 
(89),  and  (90)  is  decomposed  into 


/(2) 

uz 


p  w  +  U  w 
x  xz  My  yz 


so  that 


wxz s  ^30 +  2Q*2i  +  ^12 +  2<xP  +  yQ)^2o +  2(xQ  +  yR)*n 

+  (Px2  +  2Qxy  +  Ry2k 
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wyz  =  ^21  +  2Q^i2 +  r*o3  +  2(xP  +  yQ)<hi  +  2(xq  +  yR)^02 

+  (Px2  +  2Qxy  +  Ry2)^01 

The  ^mn  are  defined  by  equations  (96)  through  (105).  Then 

>\  n,  ... 

’(F)  -  it! 
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w  uR 


(147) 

given  by 

(148) 

(149) 

(150) 
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(150) 


9.2.2  Changes  in  On-Body  Panel  Formulas  for  Intermediate  and  Far  Fields 

In  all  cases,  Near-Field  and  Intermediate-Field  formulas  are  evaluated 
in  panel  coordinates.  For  the  vorticity  effects  discussed  in  this  section 
the  Far-Field  formulas  are  also  evaluated  in  panel  coordinates  by  truncating 
the  Intermediate-Field  expansions  with  the  first  terms.  The  ranges  of 
distance  from  the  panel  origin  where  the  various  formulas  are  used  remains 
unchanged. 


Thus  the  superscript  zero  components  are  evaluated  using  expansions  (70) 
for  the  source  velocities  (61). 

Also,  superscript  one  components  are  evaluated  using  expansions  (85), 
together  with  (87)  and 

Of  -  *4>f0  “  *20 
J2  =  y**0  ~  *fl 

(151) 

J3  =  x*0l  “  *11 
J4  =  y*0l  ~  *02 

In  an  exception  to  the  above-stated  rule,  expansion  (85)  retains  two 
terms  in  the  Far-Field  if  p  +  q  =  1.  This  is  to  ensure  that  expansions 
are  terminated  with  area  moments  of  even  order. 


For  the  superscript  two  components,  equations  (150)  are  used  with 
from  equation  (106)  with  V^R  from  (108)  and  with  (149)  replaced  by 

wxz  =  'xVo)R  +  3z(P830  +  2Qe21  +  Re12^ 

wyz  =  "yVujR  +  3z(p021  +  2Q012  +  R603^ 


mn 


(152) 
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In  equation  (106b)  both  terms  are  retained  in  the  Far-Field  for  the  reason 
stated  above. 


9.2.3  Changes  in  the  Formulas  Needed  to  Obtain  the  Effect  of  a  Wake  Panel 

The  wake  is  characterized  by  the  fact  that  the  equivalent  dipole  distribu¬ 
tion  is  constant  along  N-lines.  This  necessitates  changes  in  the  basic  assembly 
formulas  (145),  (147)  and  (150).  However,  the  formulas  for  the  individual 
quantities  appearing  in  these  equations  in  the  Near-,  Intermediate-,  and 
Far-Fields  are  unchanged. 


The  total  arc  length  of  an  N-line  from  lower-surface  trailing  edge  to 
upper-surface  trailing  edge  is  denoted  L  (total)  in  reference  1  and  given 
by  equation  (7.9.2).  This  quantity,  evaluated  on  the  first  and  second  N-lines 
of  the  panel,  is  used  in  the  wake  formulas. 


Equations  (145)  are  replaced  by 
V<;>(S)  =  0 


,(0), 

WX 


Vv„  (F)  =  0 


V^(S)  =  (145)  with  h^  replaced  by  L<j  (total) 

V^(F)  =  (145)  w.ith  hp  replaced  by  Lp  (total) 

Lc  (total) 

- c(n-|  +  n,)| 


V(°}(S)  =  - 

U)Z  ' 


V(°)( F)  = 
mz  ' 


w 

.p  (total) 


Vy  (source) 


w 


-  c(n-.  +  n-j) 


3'J  y 


V  (source) 


(153) 


Equations  (147)  are  replaced  by 

l 


vll’is)  .  2 
vO’lF)  =  -2 


s  (total) 
w 

Lp  (total) 


w 


-c(ni  +  n3)|(PJ2*  +  Q0$) 
(PJ|  +  QJ$) 


-  c(n-|  +  n3) 


^’(S)  -  2'*01  - 
-  -2<*0,  ♦ 


^-c(n-,  +  n3)  (H  +  2PJf  +  2QJ*) 


w 


c(^  +  n3)  (H  +  2PJJ  +  2QJ*) 


(154) 
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(154) 


V<’)(S)  =  -20$ 


»">  <F>  -  *2cJ4 


Equations  (150)  are  replaced  by 


V<2)(S)  =  0 


V^(F)  =  o 

b)X  '  ' 


(21 

^uy'S)  =  unchan9ed  with  h^  replaced  by  (total) 
(21 

V^y '(F)  =  unchanged  with  hp  replaced  by  Lp  (total) 


V(2)(S)  = 

<j)Z 


(total) 


w 


-  c(n-j  +  n3) 


w 

yz 


V(2)(F)  =  - 

coZ 


Lp  (total) 


-  c(n-|  +  n3) 


w 

yz 


(155) 


The  option  exists  to  make  the  last  panel  of  the  wake  semi-infinite  in 
the  manner  of  reference  1 . 


9.3  Modification  for  Parabolic  Chordwise  Vorticity  Variation 

As  discussed  in  reference  5,  the  assumption  of  constant  vorticity 
around  a  wing  section  can  lead  to  numerical  difficulties  on  certain  wings 
with  very  thin  trailing  edges.  Accordingly,  a  second  chordwise  variation 
option  is  required.  The  important  consideration  is  to  have  the  "bound" 
vorticity  strength  approach  zero  at  the  trailing  edge  on  both  upper  and  lower 
wing  surfaces.  This  is  accomplished  by  a  quadratic  global  variation  of 
vorticity  as  a  function  of  arc  length  along  an  N-line.  While  only  two  global 
chordwise  variations  have  been  incorporated  into  the  present  method,  many 
such  variations  are  possible.  If  in  the  future  numerical  difficulties  should 
be  encountered  for  certain  types  of  geometry,  say,  for  example,  a  wing  that 
is  very  thin  over  its  entire  length,  a  properly  chosen  chordwise  vorticity 
variation  could  eliminate  these  difficulties.  Moreover,  as  will  be  seen 
below,  the  required  modifications  to  the  program  are  quite  minor.  This 
flexibility  is  due  to  the  use  of  vorticity  as  an  auxiliary  singularity. 
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To  implement  the  parabolic  vorticity  option,  the  linear  variation  of  the 
equivalent  dipole  strength  along  an  N-line  (138)  is  replaced  by  a  cubic  vari¬ 
ation  having  zero  derivative  at  the  upper  and  lower  trailing  edge.  Specifically, 
(138)  is  replaced  by 


u 


Bs  3 


L  (total) 


-  2 


L  (total) 


0  56) 


where  all  symbols  have  the  same  meaning  as  in  section  9.2.  The  above  is  a 
global  variation.  The  variation  over  an  individual  panel  can  be  no  higher  a 
degree  than  quadratic,  and  in  the  initial  version  of  the  present  method  has 
been  taken  as  linear.  It  is  assumed  that  the  dipole  distribution  on  a  panel 
agrees  with  (156)  at  the  corners  of  the  panel  and  varies  linearly  in  between. 
Thus,  the  overall  behavior  is  that  of  an  inscribed-polygon  approximation  to 
(156). 


For  an  individual  panel  the  arc  length  measured  along  an  N-line  from  the 
trailing  edge  to  the  lower  corner  of  the  panel  is  L  (equation  (7.2.23)  of 
reference  1),  while  the  arc  length  associated  with  the  upper  corner  is 
L  +  d  (equation  (7.2.22)  of  reference  1).  Both  L  and  d  are  subscripted 
with  F  or  S  to  denote,  respectively,  the  first  or  second  N-lines.  The 
linear  function  that  agrees  with  (145)  at  these  two  values  of  arc  length  is 

v  -  B(H  +  10  (157) 

where  on  the  two  N-lines  the  constants  H  and  I  have  the  values 
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where  all  symbols  have  the  same  meaning  as  in  section  9.2.  Thus  the  quadratic 
form  (140)  for  the  variation  of  dipole  strength  over  a  panel  is  replaced  by 


_  Vf  Vs  _  .  BFHF  BSHS  .  BSISnl  BFIFn3  .  .  BsVl  ~  BFHFn3 

y - z - Cn  + - - - n  + - - - K  + - 


w 


w 


w 


w 


+  c(Bp  -  Bs)(n  -  n3)(n  -  tvj)  +  d(Bp  - 


(159) 


where  in  the  initial  version  of  the  present  method  c  and  d  are  set  equal 
to  zero.  The  dipole  derivative  formulas  (141)  are  modified  in  an  obvious  way. 


Thus  the  elaborate  vorticity  velocity  formulas  of  section  9.2  are  modi¬ 
fied  to  account  for  parabolic  vorticity  variation  by  the  following  procedure: 

a.  Terms  containing  c  or  d  are  not  changed. 

b.  In  terms  containing  hp  or  h$  these  quantities  are  replaced  by 
Hp  or  H<> ,  respectively 

c.  All  other  terms  are  multiplied  by  Ip  or  1^  as  appropriate. 

In  the  constant  chordwise  vorticity  option,  the  parameter  c  is  nonzero  on 
wake  panels  if  the  "piecewise  linear"  spanwise  variation  of  vorticity  is 
used  (section  7.9  of  reference  1).  However,  if  the  parabolic  chordwise 
vorticity  option  is  used,  c  is  taken  as  zero  on  all  panels. 


9.4  The  Piecewise  Linear  Spanwise  Vorticity  Variation 

The  third  aspect  of  the  global  vorticity  algorithm  is  the  variation 
from  N-line  to  N-line  across  the  "span"  of  a  panel.  The  base  method  of 
reference  1  has  two  options:  constant  and  linear.  In  the  former  the  equiv¬ 
alent  dipole  strength  is  discontinuous  across  an  N-line  by  an  amount  equal  to 
the  product  of  the  spanwise  gradient  and  the  panel  width.  Concentrated 
trailing  vortex  filaments  of  this  strength  lie  along  the  N-line.  This  repre¬ 
sentation  is  inconsistent  with  a  higher-order  formulation  and  has  been 
excluded  from  the  vorticity  velocity  formulas  of  section  6.3.  Accordingly, 
the  present  method  has  only  one  type  of  spanwise  variation,  the  linear  one. 

Once  the  vorticity  velocities  V^F)  and  V  (R)  have  been  calculated 
by  the  procedures  of  sections  9.1,  9.2  and  9.3,  they  may  be  used  in  exactly 
the  same  way  as  the  corresponding  quantities  of  the  base  method.  Thus  the 
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procedure  of  section  7.11  of  reference  1  may  be  used  without  change  to 
generate  the  vorticity  onset  flows  that  are  used  to  satisfy  the  Kutta 
condition. 


1 


10.0  THE  KUTTA  CONDITION 

The  Kutta  condition  is  applied  in  one  of  two  ways:  a  condition  of  flow 
tangency  in  the  wake  or  a  condition  of  pressure  equality  on  the  upper  and 
lower  surfaces  of  the  trailing  edge,  which  amounts  to  a  condition  of  equal 
velocity  magnitude,  (section  7.13,  reference  1).  The  wake  tangency  condition 
is  applied  in  exactly  the  same  manner  as  reference  1.  This  section  is  con¬ 
cerned  with  the  equal-pressure  condition,  which  uses  a  different  iterative 
procedure  from  that  of  reference  1 . 

As  a  numerical  approximation,  the  Kutta  condition  may  be  applied  by 
equating  pressures  at  the  control  points  of  the  two  panels  adjacent  to  the 
trailing  edge  on  the  upper  and  lower  surfaces  of  the  wing.  Alternatively, 
velocities  at  the  upper  surface  control  points  of  the  few  panels  nearest  the 
trailing  edge  could  be  extrapolated  to  obtain  velocity  components  "at"  the 
trailing  edge  upper  surface,  and  the  same  could  be  done  for  the  lower  surface. 
This  last  would  allow  application  of  the  Kutta  condition  more  nearly  at  the 
trailing  edge,  and  the  analogy  of  this  procedure  yields  considerable  improve¬ 
ment  in  accuracy  in  two-dimensional  cases.  However,  in  the  initial  version 
of  the  present  method,  pressure  equality  is  applied  at  the  control  points  of 
the  adjacent  panels. 

However  the  point  of  application  of  the  Kutta  condition  is  chosen,  the 

logic  of  the  calculation  is  the  same.  In  particular,  a  velocity  vector  can 

be  calculated  at  the  upper  and  at  the  lower  trailing  edge  point  for  each 

strip  of  panels  (figure  1)  and  for  each  onset  flow  —  both  the  uniform  stream 

and  the  several  vorticity  onset  flows  (section  9.4).  It  is  necessary  to 

introduce  some  notation.  A  subscript  m  denotes  quantities  associated  with 

the  m-th  lifting  strip  of  panels,  and  m  =  1,  2,  ...»  L,  where  L  is  the 

total  number  of  lifting  strips.  Superscript  ®  denotes  velocities  obtained 

in  the  flow  solution  for  a  uniform  onset  flow,  and  superscript  k  =  1,  2,  ...»  L 

denotes  velocities  obtained  in  the  flow  solution  due  to  the  vorticity  onset 

ikl 

flow  associated  with  the  K-th  lifting  strip.  Bv  '  denotes  the  unknown 
multiplicative  constant  to  be  applied  to  the  vorticity  solution  correspond¬ 
ing  to  the  k-th  lifting  strip.  Finally,  (Up)  or  (Lo)  denote  velocities 
at  the  upper  or  lower  trailing-edge  point,  respectively.  Thus  the  velocity 
vectors  at  the  trailing-edge  points  of  the  m-th  strip  are 
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*„«lp>  =  'i’W)  B(k)VIJ1k>(Up) 

k=l 

m  =  1,  2,  ....  L  (160) 

Vm(Lo)  =  vW(Lo)  +2B(k)^k)(Lo) 

k=l 

These  are  analogous  to  the  velocity  vectors  (7.13.3)  of  reference  1.  A  unit 
vector  for  each  strip  of  panels  is  defined  by 

VJUP)  *  v  (Lo) 


Vm(Avg)  = 


m' 


m' 


|Vm(Up)  +  Vm(Lo)| 


m  =  1 ,  2,  ...»  L  (161 ) 


The  condition  that  the  velocity  vectors  at  the  upper  and  lower  trail ing- 
edge  points  of  each  strip  have  equal  magnitudes  may  be  written 


Vm(Up)  •  Vm(Avg)  =  Vm(Lo)  •  VjAvg) 


m  =  1,  2 . L  (162) 


Equations  (162)  consist  of  L  simultaneous  quadratic  equations  for  the  L 
fk) 

unknowns  Bv  ‘ .  All  coefficients  of  equations  (162)  can  be  expressed  in  terms 
of  the  2m(m  +  1)  individual  velocities  appearing  in  (160).  Equations  (162) 
are  solved  iteratively  beginning  with  =  0.  At  any  stage  of  the  iteration 

the  Vm(Avg)  are  evaluated  from  the  currently  known  set  of  B^.  These  are 

used  in  (162)  to  obtain  an  L  x  L  set  of  linear  equations  for  the  successive 

(kl  - 

set  of  Bx  ' .  The  iteration  has  been  thoroughly  tested  in  the  base  method  of 

reference  1  and  has  converged  very  rapidly  in  all  cases. 
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11.0  CALCULATED  RESULTS 


To  illustrate  the  improvements  that  the  higher-order  method  of  this 
report  has  been  able  to  effect  compared  to  the  base  method  of  references  1 
and  2,  a  number  of  cases  have  been  run  for  which  high-accuracy  solutions  are 
available.  While  it  would  undoubtedly  be  useful  to  run  more  such  cases,  the 
results  presented  below  strongly  indicate  the  effectiveness  of  the  present 
method,  even  in  its  initial  form. 

The  triaxial  ellipsoid  is  the  only  truly  three-dimensional  body  for 
which  an  analytic  solution  exists.  A  geometry  was  selected  having  semiaxes 
equal  to  1,  2  and  H  in  the  x,  y,  and  2  directions,  respectively,  as  shown  in 
the  sketch  of  figure  6.  Calculations  were  performed  for  two  uniform  onset 
flows  -  one  parallel  to  the  x-axis  and  one  to  the  z-axis  -  using  16  strips 
of  panels  defined  by  section  cuts  (N-lines)  at  constant  values  of  y  from 
+2  to  -2.  Each  strip  had  30  panels  around  the  elliptical  cross-section  — 
a  total  of  480  effective  panels  (symmetry  was  employed).  Calculated  and 
analytic  surface  velocities  along  the  curve  in  the  xz-plane  are  compared  in 
figure  6.  As  pointed  out  in  references  5,  6  and  7,  the  base  method  enjoys 
a  certain  amount  of  cancellation  of  errors  when  calculating  flow  about  smooth 
convex  shapes.  Thus  the  higher-order  procedure  cannot  be  expected  to  yield 
consistent  gains  in  accuracy  for  such  cases.  In  view  of  this,  the  small  but 
definite  improvement  in  accuracy  exhibited  by  the  higher-order  results  in 
figure  6  is  quite  significant. 

Greater  gains  in  calculational  accuracy  can  be  realized  by  the  higher- 
order  method  when  the  body  has  concave  regions.  The  geometry  selected  is  the 
circular  nacelle  shown  in  the  sketch  of  figure  7.  With  the  onset  flow  parallel 
to  the  axis  of  symmetry,  the  total  flow  is  axisymmetric  and  a  high-accuracy, 
graphically  exact  solution  can  be  obtained  by  the  method  of  reference  6. 
Nevertheless,  this  is  a  meaningful  case  for  comparing  the  three-dimensional 
base  and  higher-order  methods,  because  they  do  not  take  advantage  of  the 
axial  symmetry  but  panel  the  nacelle  both  axially  and  circumferentially.  The 
three-dimensional  results  shown  in  figure  7  were  obtained  using  a  distribution 
of  17  panels  axially  on  the  inner  surface  of  the  nacelle  and  the  same  number 
on  the  upper  surface.  Circumferentially,  18  effective  panels  were  employed, 
so  that  each  panel  subtends  an  angle  of  20°.  This  is  a  lifting  case.  Kutta 
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conditions  are  applied  along  the  trailing  edge  of  the  nacelle,  and  the 
vorticities  on  the  strips  are  adjusted  accordingly.  Because  the  axi symmetric 
solution  of  reference  6  employed  a  much  larger  number  of  panels  axially  than 
either  of  the  three-dimensional  methods,  it  applied  the  Kutta  condition  of 
equal  upper  and  lower  surface  pressures  at  a  point  nearer  the  true  trailing 
edge.  This  disturbs  the  agreement  in  the  trail ing-edge  region,  but  has 
relatively  little  effect  at  forward  locations.  Two-dimensional  experience 
indicates  that  this  situation  can  be  alleviated  either  by  extrapolating  the 
Kutta  condition  to  the  trailing  edge  or  by  using  small  panels  adjacent  to  the 
trailing  edge.  Figure  7  dramatically  illustrates  the  effect  of  concavity  on 
the  exterior  surface,  where  the  body  is  convex,  all  three  calculated  solutions 
agree  quite  well  (away  from  the  trailing  edge).  On  the  concave  interior  surf¬ 
ace  the  lower  order  method  of  reference  1  gives  results  that  are  seriously  in 
error,  while  the  higher-order  method  of  the  present  report  agrees  with  the 
high-accuracy  solution  almost  exactly. 

The  most  impressive  improvements  in  accuracy  due  to  inclusion  of  higher- 
order  effects  occur  in  interior  flows.  To  illustrate  this,  the  geometry 
selected  is  the  circular  duct  of  area  ratio  four  shown  in  figure  8.  The  duct 
was  paneled  using  strips  of  panels  around  the  circular  cross  sections.  Each 
strip  of  panels  has  a  width  equal  to  one- third  of  the  maximum  duct  radius  — 
a  total  of  42  strips  along  the  length  of  the  duct.  Each  strip  contains  18 
effective  panels  (symmetry  is  employed)  around  the  circular  cross  sections, 
so  that  each  panel  subtends  an  angle  of  20°.  For  purposes  of  comparison, 
a  high  accuracy,  graphically  exact  solution  for  the  duct  has  been  obtained  by 
using  the  higher-order  axisymmetric  method  of  reference  6.  Several  calcula¬ 
tions  have  been  made,  each  with  twice  as  many  panels  along  the  duct  as  the 
previous  calculation.  This  study  has  verified  that  the  solution  obtained  by 
the  higher-order  axisymmetric  program  using  140  panels  along  the  duct  may  be 
considered  graphically  exact,  and  accordingly  it  is  used  as  a  standard  of 
comparison  in  figure  8.  The  results  of  three  different  three-dimensional 
calculations  are  compared  with  the  high  accuracy  solution  in  figure  8.  The 
first  two  use  the  42  x  18  panel  distribution  described  above  -  a  total  of 
756  effective  panels.  One  calculation  was  performed  with  the  higher-order 
method  of  the  present  report  and  one  with  the  base  method  of  references  1 
and  2.  The  base-method  solution  Is  so  Inaccurate  as  to  be  unusable,  while 
the  higher-order  solution  is  scarcely  distinguishable  from  the  exact.  To 
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provide  another  standard  of  judgment,  the  results  of  two  base-method  calcu¬ 
lations  using  a  greatly  increased  panel  number  are  also  shown  in  figure  8. 
These  cases,  which  have  been  taken  from  reference  2,  each  use  4200  effective 
panels.  One  has  42  strips  of  panels  along  the  duct,  each  with  100  panels 
around  the  circular  cross  section.  The  other  has  84  strips  of  panels  along 
the  duct,  each  with  50  panels  around  the  circular  cross  section.  The  calcu¬ 
lated  surface  velocities  for  both  cases  agree  fairly  well  (reference  2), 
and  so  they  are  shown  as  a  single  curve  in  figure  8.  Even  with  this  greatly 
increased  panel  number,  the  errors  in  the  base-method  calculation  are  an  order 
of  magnitude  larger  than  those  of  the  higher-order  calculation.  Since  approx¬ 
imately  5.5  times  as  many  panels  are  employed,  the  computing  time  for  the 
large  panel  number  case  of  figure  8  is  at  least  20  times  larger  than  that  of 
the  higher-order  method.  Thus,  for  the  duct  of  figure  8,  the  higher-order 
method  obtained  a  factor  10  increase  in  accuracy  and  a  factor  20  decrease  in 
computing  time. 


12.0  PRINCIPAL  NOTATION 


A,B  coefficients  of  linear  terms  in  the  expression  for  local  surface 
shape,  equation  (110) 

n  unit  normal  vector  with  components  n„,  n  ,  n, 

C  n  C 

P,Q,R  coefficients  of  quadratic  terms  In  the  expression  for  local  surface 
shape,  equations  (6)  or  (110).  These  "local  curvatures"  are  funda¬ 
mental  to  the  higher-order  procedure. 

r  distance  between  a  field  point  (x,y,z)  and  a  point  U,n,t)  on 

the  curved  panel 

rf  distance  between  a  field  point  (x,y,z)  and  a  point  U»n,0)  on 

the  projected  flat  panel 

rQ  distance  from  a  field  point  (x,y,z)  to  the  control  point  of  a 

curved  panel 

t  panel  dimension.  Usually  the  maximum  chord  of  the  projected  flat 

panel . 

V  fluid  velocity.  Used  with  various  subscripts  and  superscripts  to 

denote  components  and  velocities  due  to  various  effects. 

x,y,z  Cartesian  coordinates  of  a  point  in  space,  usually  in  the  panel 

coordinate  system.  Used  as  subscripts  to  denote  partial  derivatives 
and  components  of  vectors. 

the  summed  quadratic  terms  of  the  local  shape,  equation  (6).  Sub¬ 
scripts  denote  derivatives. 

v  dipole  strength.  Used  with  various  subscripts  to  denote  derivatives 

at  the  origin  of  panel  coordinates 

C.n.c  Cartesian  coordinates  of  a  point  on  the  curved  panel  in  its  own 

coordinate  system.  5  and  n  denote  coordinates  of  a  point  on  the 
projected  flat  panel.  Used  as  subscripts  to  denote  derivatives 
and  vector  components. 

source  density.  Used  with  various  subscripts  to  denote  derivatives 
at  the  origin  of  panel  coordinates. 
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velocity  potential.  Used  with  various  subscripts  and  superscripts 
to  denote  potentials  due  to  various  effects. 

vorticity  strength 

The  equations  of  this  report  are  not  intended  to  stand  Independent  of 

those  of  reference  1.  Similarly,  quantities  defined  there  are  not 

repeated  here.  These  include:  I.  *  ,  J  ,  u  (with  various 

mn  mn  mn 

subscripts),  m32,  b32,  B,  L,  w,  h,  a^,  sK>  nK  and  subscripts  F 
and  S. 
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Figure  3.  A  general  curved  surface  panel  and  its  projection  in  the  tangent 
plane. 


STRIP  CONTAINING 


Figure  4.  The  input  points  of  adjacent  panels  used  to  compute  local 
surface  curvatures. 
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PANEL  IN  QUESTION; 

iiitij) . 


4 


Figure  5.  Adjacent  panels  used  in  the  least-squares  fit  for  source  density 
derivatives. 


Figure  6.  Comparison  of  calculated  and  analytic  surface  velocities  on  a 
triaxial  ellipsoid  with  axes  ratios  1,  2,  1/2  in  the  x,y,z 
directions.  Velocities  shown  are  along  the  curve  lying  in  the 
xz-plane. 
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Figure  7.  Calculated  surface  pressure  distribution  on  a  circular  nacelle. 


Calculated  surface  velocity  distributions  on  a  circular  duct  of  area  ratio  four 


